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ABSTRACT 

Robarge,  Tyler  W.  MSAE,  Purdue  University,  August,  2005.  Laminar  Boundary- 
Layer  Instabilities  on  Hypersonic  Cones:  Computations  for  Benchmark  Experiments. 
Major  Professor:  Steven  P.  Schneider. 

Although  significant  advances  have  been  made  in  hypersonic  boundary-layer  tran¬ 
sition  prediction  in  the  last  several  decades,  most  design  work  still  relies  on  unreli¬ 
able  empirical  correlations  or  wind-tunnel  tests.  Codes  using  the  semi-empirical  eN 
method  will  need  to  be  well  verified  and  validated  before  being  used  for  expensive 
flight  vehicles.  The  code  package  STABL  and  its  PSE-Chem  stability  solver  are 
used  to  compute  first  and  second-mode  instabilities  for  both  sharp  and  blunt  cones 
at  wind-tunnel  conditions  using  a  Navier-Stokes  mean-flow  solution.  Computations 
are  performed  for  Stetson’s  3.81  mm  nose-radius  cone,  a  sharp  cone  at  Mach  3.5, 
a  large-bluntness  cone  at  Mach  8,  and  sharp  and  blunt  cones  corresponding  to  the 
experiments  of  Rufer.  Comparisons  to  previous  computations  by  other  researchers 
show  differences  on  the  order  of  10%  in  local  amplification  rates  and  frequencies, 
but  better  agreement  is  obtained  for  the  transition  location.  Many  issues  are  ex¬ 
amined  for  verification  and  validation,  including  the  laminar  transport  properties, 
the  freestream  boundary  conditions,  and  the  effect  of  freestream  thermal  nonequi¬ 
librium.  This  work  helps  to  verify  and  validate  STABL,  extend  its  applicability  to 
low-temperature  flows,  and  develop  the  methodologies  for  using  STABL. 
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1.  Introduction 

Despite  more  than  fifty  years  of  hypersonics  research,  many  critical  areas  are  still 
poorly  understood  [1].  Hypersonic  flight  remains  a  top  priority  in  defense  research, 
as  any  vehicle  that  ascends  to  or  from  space  must  traverse  the  hypersonic  regime.  An 
Air  Force  Scientific  Advisory  Board  study  found  that  hypersonics  will  be  necessary 
to  fulfill  the  Air  Force’s  vision  of  “controlling  and  exploiting  the  full  aerospace  con¬ 
tinuum”  [2],  The  reduction  in  forward-deployed  forces  since  the  end  of  the  Cold  War 
has  increased  the  need  for  hypersonics,  as  a  need  remains  to  be  able  to  strike  high- 
value,  time-critical  targets  anywhere  in  the  world  in  minutes  or  hours  [3].  Hypersonic 
flight  is  a  key  component  of  most  plans  to  provide  that  capability. 

The  location  and  extent  of  laminar-to-turbulent  boundary-layer  transition  is  a 
critical  parameter  in  hypersonic  vehicle  design  [1],  The  transition  location  directly 
affects  estimates  of  aeroheating  and  skin  friction  drag,  which  in  turn  affect  heat 
shield  weight  and  materials,  range,  and  payload  capacity.  Infrared  and  radar  sig¬ 
natures,  as  well  as  communications  with  the  vehicle,  are  also  affected  by  transition. 
Vehicle  designers  can  now  compute  the  laminar  and  turbulent  aeroheating  with  good 
accuracy;  in  many  cases,  the  largest  uncertainty  in  calculating  the  total  heat  flux  to 
the  vehicle  results  from  the  estimates  of  transition  location  [4],  As  an  example,  the 
Defense  Science  Board  found  that  estimates  for  transition  location  on  the  National 
Aerospace  Plane  ranged  from  20%  to  80%  of  the  body,  and  the  estimate  used  could 
affect  the  vehicle  gross  takeoff  weight  by  a  factor  of  two  or  more  [5].  Transition  is  es¬ 
pecially  important  on  vehicles  that  spend  a  significant  amount  of  time  in  high-speed, 
high-altitude  atmospheric  flight,  such  as  hypersonic  cruise  vehicles  and  manuevering 
reentry  vehicles. 

Figure  1.1  shows  a  magnified  portion  of  a  shadowgraph  of  a  sharp  cone  model  in 
free  flight  at  Mach  4.3  in  the  Naval  Ordinance  Lab  ballistics  range.  This  photograph, 
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Figure  1.1:  Shadowgraph  of  a  sharp  cone  in  free  flight  at  Mach  4.3.  Picture  first 
published  by  Schneider  [4], 


first  published  by  Schneider  [4],  illustrates  many  of  the  salient  features  of  hypersonic 
boundary-layer  transition.  The  boundary  layer  on  the  upper  surface  of  the  cone  is 
mostly  laminar  with  discrete  turbulent  spots.  Small  shocks  are  visible  forward  of  the 
spots.  The  visible  lower  surface  is  entirely  turbulent,  and  acoustic  noise  can  be  seen 
radiated  from  the  boundary  layer.  A  vast  difference  in  the  boundary-layer  thickness 
between  the  laminar  and  turbulent  states  is  evident. 

1.1  Transition  Prediction  Methods 

The  overall  transition  process  is  not  well  understood  [6].  Environmental  distur¬ 
bances,  which  could  come  from  either  the  freestream  or  the  body,  enter  the  boundary 
layer  through  a  process  called  receptivity.  These  disturbances  amplify  through  one 
or  more  of  many  possible  mechanisms.  This  thesis  will  consider  only  first  and  sec¬ 
ond  mode  amplification,  but  roughness-induced  transition,  cross-flow  vortices,  and 
Gortler  instabilities  are  all  examples  of  other  mechanisms  that  are  dominant  under 
various  conditions.  As  the  disturbances  amplify,  they  will  eventually  breakdown  in 
a  nonlinear  fashion,  form  turbulent  spots,  and  develop  into  full  turbulence. 
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Despite  this  complex  nature  of  transition,  transition  prediction  for  vehicle  design 
has  historically  relied  on  empirical  correlations  and  extrapolations  from  wind  tunnel 
experiments.  A  large  body  of  empirical  and  semi-empirical  correlations  exist;  many 
of  these  were  cataloged  and  evaluated  by  Berkowitz  et  al.  [7].  One  popular  correlation 
is  shown  in  Equation  1.1. 


Reg  _ 


(1.1) 


The  constant  C  is  specified  based  on  empirical  data  and  is  usually  in  the  range  of 
100-500.  Transition  is  assumed  to  occur  when  the  local  boundary-layer  properties 
satisfy  Equation  1.1.  The  designers  of  the  Space  Shuttle  used  Equation  1.1  and 
applied  a  correction  based  on  the  surface  roughness  caused  by  misaligned  Thermal 
Protection  System  tiles  [8]. 

As  Berkowitz  et  al.  showed,  many  correlations  exist,  and  most  fit  some  subset 
of  the  experimental  data.  However,  no  empirical  model  can  accurately  predict  tran¬ 
sition  for  a  general  dataset  [1].  In  fact,  Schneider  [9]  has  shown  that  for  a  general 
data  set  the  predictions  generated  by  different  correlations  scatter  by  a  factor  of  3 
in  momentum  thickness  Reynolds  number  and  an  order  of  magnitude  in  arc  length 
Reynolds  number.  None  of  the  correlations  can  be  expected  to  provide  accurate 
predictions  for  configurations  and  conditions  outside  of  the  database  used  for  their 
development. 

Extrapolating  from  wind  tunnel  experiments  is  also  problematic.  Transition  de¬ 
pends  on  a  multitude  of  factors,  including  local  Mach  number,  flow  enthalpy,  unit 
Reynolds  number,  surface  roughness,  reacting  chemistry,  and  mass  ablation  [6].  Any 
single  ground  facility  can  only  match  a  small  subset  of  these  parameters,  and  even 
then  determining  the  location  and  extent  of  transition  in  an  accurate  and  consistent 
manner  is  not  a  trivial  task  [1],  Further  complicating  matters,  the  acoustic  noise 
produced  by  the  turbulent  boundary  layer  on  the  walls  of  conventional  hypersonic 
wind  tunnels  produces  a  freestream  disturbance  environment  markedly  different  from 
that  of  flight.  These  disturbances  have  been  shown  to  affect  not  only  the  transition 
location,  but  also  the  parametric  trends  [10].  This  has  led  to  a  push  for  the  devel- 
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opment  of  ‘quiet’  tunnels,  but  that  itself  is  an  exercise  in  boundary-layer  transition 
prediction  and  control.  There  are  currently  no  hypersonic  wind  tunnels  that  are 
quiet  at  significant  Reynolds  numbers,  anywhere  in  the  world  [11]. 

Given  the  constraints  of  ground  testing,  the  only  situation  that  can  accurately 
capture  all  of  the  relevant  flow-field  characteristics  is  flight.  However,  flight  testing 
usually  only  permits  surface  measurements,  and  recovery  of  the  actual  vehicle  is 
often  impossible.  Flight  testing  is  extremely  expensive,  and  sufficient  funds  are  not 
likely  to  be  available  in  the  forseeable  future. 

Much  research  in  hypersonic  bound  ary- layer  transition  has  been  accomplished  in 
the  last  fifty  years.  In  particular,  Linear  Stability  Theory  (LST)  [12]  and  the  Parab¬ 
olized  Stability  Equations  (PSE)  [13]  have  been  coupled  with  the  semi-empirical  eN 
method  to  predict  transition  on  simple  geometries,  such  as  cones  and  flat  plates. 
More  details  on  all  of  these  methods  will  be  provided  in  Chapter  2.  Recently,  Ma¬ 
lik  [14]  showed  that  transition  locations  could  be  correlated  on  two  hypersonic  flight 
tests  using  the  PSE  for  chemically  reacting  flows  and  the  eN  method.  Although 
these  methods  have  shown  success  for  many  applications  in  research  settings  [6],  hy¬ 
personic  vehicle  designers  are  hesitant  to  stake  expensive  programs  on  these  newer 
methods  until  they  have  been  shown  to  be  better  validated,  faster,  and  simpler 
to  use.  Given  these  constraints,  Oberkampf  and  Blottner  [15]  assert  that  reliable 
engineering  techniques  for  transition  prediction  are  currently  lacking. 

Since  hypersonic  vehicles  remain  a  priority,  and  current  design  methods  are  in¬ 
adequate  for  many  applications,  an  improved  method  for  hypersonic  boundary-layer 
transition  prediction  is  needed.  The  overall  project  that  this  work  is  a  part  of  aims 
to  provide  mechanism-based  methods  that  are  suitable  for  design  purposes,  includ¬ 
ing  predictions  using  the  PSE  and  the  eN  method.  The  present  work  is  concerned 
primarily  with  the  verification  and  validation  of  the  2D/axisymmetric  version  of  the 
PSE/e^  code  package.  This  code  package  will  not  be  the  final  solution  to  the  prob¬ 
lem  of  transition  prediction,  as  the  only  mechanisms  currently  included  are  the  first 
and  second  mode,  and  the  effects  of  receptivity  and  nonlinear  breakdown  are  treated 
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empirically.  However,  it  is  hoped  that  the  inclusion  of  more  of  the  relevant  flow 
physics  will  result  in  a  more  accurate  and  more  robust  method  than  those  which  are 
currently  in  use. 

1.2  Verification  and  Validation 

Before  computational  results  can  be  used  to  make  important  decisions,  measures 
must  be  taken  to  ensure  confidence  in  the  accuracy  of  the  predictions.  Although  this 
issue  has  received  considerable  attention  in  recent  years,  developments  in  this  area 
have  not  kept  pace  with  increases  in  code  size  and  complexity  and  the  increasing 
reliance  on  computations  to  reduce  the  number  and  scope  of  wind  tunnel  and  flight 
tests  [16].  Attention  to  code  quality,  implementation,  and  the  applicability  of  the 
algorithms  and  assumptions  is  required  throughout  all  stages  of  code  creation  and 
use,  but  even  this  does  not  ensure  that  a  code  is  free  from  defects.  As  an  example, 
an  inexperienced  user  can  get  incorrect  results  from  a  validated  code. 

Blottner  [17]  defines  code  validation  as  solving  the  right  governing  equations  and 
code  verification  as  solving  the  governing  equations  right.  Examples  of  validation 
include  confirming  the  accuracy  of  the  assumptions  made,  such  as  when  using  the 
Euler  instead  of  the  Navier-Stokes  equations,  or  comparing  the  accuracy  of  the  clo¬ 
sure  models  throughout  the  relevant  parameter  space.  Verification  means  evaluating 
the  accuracy  of  the  numerical  procedures  used  to  solve  the  governing  equations.  Ex¬ 
amples  of  verification  include  checking  for  general  programming  errors  and  ensuring 
that  results  are  converged  and  grid  independent.  Code  validation  can  only  be  ac¬ 
complished  by  detailed  comparisons  to  experimental  data  [18],  and  the  validation  of 
applicable  equations  and  models  is  a  responsibility  of  the  community  at  large.  In 
general,  a  code  must  be  verified  before  it  can  be  fully  validated. 

Oberkampf  and  Blottner  [15]  divide  the  sources  of  physical  modeling  errors  into 
three  categories:  the  partial  differential  equations  (PDEs),  the  closure  models,  and 
the  boundary  conditions.  The  set  of  PDEs  analyzed  must  be  appropriate  for  the 
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flow  phenomena  of  interest.  In  addition,  the  accuracy  of  the  solution  depends  not 
just  on  the  continuum  mathematical  model  used  to  describe  the  flow,  but  also  on 
the  accuracy  of  the  discrete  approximation  used  to  solve  that  model.  The  closure 
models  must  be  accurate  and  valid  for  the  full  parameter  space  of  temperature,  den¬ 
sity,  etc.  The  boundary  conditions  are  arguably  the  hardest  to  accurately  model,  as 
this  requires  detailed  knowledge  of  the  inflow  and  outflow  profiles  and  disturbance 
spectra.  This  is  rarely  available,  even  in  careful  experiments  designed  for  code  val¬ 
idation.  Aeschliman  and  Oberkampf  [19]  provide  a  good  description  of  the  type  of 
experimental  and  computational  procedure  necessary  to  rigorously  validate  a  code. 

Verification  is  a  difficult  task,  and  it  is  never  fully  completed.  Hatton  and 
Roberts  [20]  found  from  a  study  of  seismic  data-processing  software  packages  that 
numerical  disagreement  among  commercially  available  codes  grows  at  about  1%  in 
average  absolute  difference  per  4000  lines  of  code,  and  that  the  Fortran  codes  had  a 
static  fault  rate  of  6  per  1000  lines  of  code.  These  errors  among  the  various  packages 
examined  led  to  results  that  were  all  reasonable  but  subtly  different,  making  it  im¬ 
possible  for  even  trained  scientists  to  evaluate  which  results  were  the  most  accurate. 
Given  that  these  results  were  obtained  for  a  mature,  highly  competitive  field,  where 
the  algorithms  are  well  estabilished  and  over  fifteen  independent  packages  exist  to 
perform  the  same  function,  it  is  highly  probable  that  the  situation  is  worse  for  hyper¬ 
sonic  transition  prediction,  where  funding  is  insufficient  to  support  many  competing 
products. 

Several  approaches  exist  for  code  verification,  including  grid  refinement  studies, 
varying  the  boundary  conditions,  solving  problems  that  have  an  analytical  solution, 
and  comparing  results  with  previously  verified  codes  [16,18].  Grid  refinement  studies 
and  comparisons  with  other  codes  will  be  used  in  the  present  effort.  This  work 
expands  on  previous  verification  and  validation  work  by  Johnson  and  Johnson  et 
al.  [21-24], 
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2.  Code  Description 

The  computations  presented  in  this  thesis  were  performed  using  the  STABL  code 
package  [24],  STABL  is  designed  to  be  a  comprehensive  boundary-layer  stability- 
analysis  and  transition-prediction  tool,  and  it  includes  a  mean  flow  code,  a  stability 
code,  and  other  utilities  such  as  a  Graphical  User  Interface  (GUI)  and  a  grid  gen¬ 
erator.  STABL  is  modular  in  nature,  allowing  individual  components  to  be  used 
instead  of  the  whole  package. 

2.1  History 

The  heart  of  STABL  is  its  stability  code,  PSE-Chem.  PSE-Chem  was  written 
by  Heath  B.  Johnson  as  part  of  his  Ph.D.  research  under  Professor  Graham  Candler 
at  the  University  of  Minnesota  in  the  late  1990s.  The  original  PSE-Chem  was  a 
research  code  with  a  very  utilitarian  user  interface  and  very  little  documentation 
apart  from  the  extensive  comments  within  the  source  code.  The  Navier-Stokes  code 
DPLR2D  [25],  written  by  Michael  Wright,  another  of  Professor  Candler’s  former 
students,  was  used  for  the  mean  flow  calculations. 

After  completing  his  degree,  Johnson  began  improving  the  PSE-Chem  code  and 
the  other  components  to  add  new  capabilities  and  to  make  them  more  user-friendly. 
When  the  author  began  using  the  code  package  in  December  of  2003,  he  was  the 
first  individual  outside  of  Professor  Candler’s  research  group  to  use  the  code.  At  the 
time,  the  user’s  manual  was  approximately  fifteen  pages  long,  a  minimal  GUI  was  in 
place,  and  compiling  the  codes  required  manually  editing  the  makefiles  to  configure 
system-specific  variables. 

Over  the  course  of  the  next  eighteen  months,  the  author  served  as  the  primary 
tester  for  new  components  and  features  released  by  Dr.  Johnson.  Many  improve- 
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ments  were  made  to  the  package,  allowing  for  a  much  easier  learning  curve  for  new 
users.  A  large  percentage  of  these  improvements  directly  resulted  from  suggestions 
or  ideas  raised  by  the  author.  Over  the  successive  months,  a  script  was  added  that 
allows  a  user  to  answer  a  series  of  questions  to  configure  the  system-specific  vari¬ 
ables,  which  mostly  eliminated  the  need  to  edit  the  obscure  makefiles.  The  GUI  was 
expanded  and  improved,  such  as  through  the  addition  of  a  flow  analysis  tool  to  allow 
for  easier  interpretation  of  mean  flow  results.  Algorithm  changes,  detailed  in  later 
chapters,  were  made  to  provide  more  accurate  and  robust  results. 

In  March  of  2004,  version  1.1  was  released,  and  in  April  of  2004,  the  package 
was  renamed  to  Stability  and  Transition  Analysis  for  hypersonic  Boundary  Layers 
(STABL).  The  computations  presented  in  this  thesis  were  performed  using  versions 
of  STABL  ranging  from  1.26  to  1.29.4.  The  ongoing  development  of  STABL  made 
obtaining  all  results  with  a  consistent  version  impractical.  Additional  details  about 
STABL  can  be  found  in  the  STABL  Program  Reference  [26]  and  in  Johnson  and 
Candler  [24], 

2.2  Major  Components  of  STABL 
2.2.1  Grid  Generator 

The  grids  used  for  the  mean  flow  solutions  were  created  using  the  grid  generator 
included  with  STABL.  The  model  surface  and  outer  boundaries  were  each  specified 
using  a  single  arc  section  and  one  or  more  line  segments.  In  the  case  of  a  sharp 
cone,  only  line  segments  were  used.  As  the  mean  flow  solver  uses  a  shock  capturing 
scheme,  the  domain  outer  boundary  had  to  be  carefully  defined  to  closely  follow  the 
shock  shape  to  minimize  the  number  of  grid  lines  crossing  the  shock.  Since  the  shock 
shape  is  not  known  a  priori ,  an  iterative  approach  was  used  in  which  a  coarse,  poorly 
fitted  grid  was  used  to  obtain  a  solution,  and  the  grid  parameters  were  adjusted  until 
the  outer  boundary  of  the  new  grid  matched  the  shock  shape  in  the  original  solution 
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with  a  small  offset.  A  second  mean  flow  solution  was  then  obtained  using  the  new 
grid  parameters  with  a  higher  mesh  density. 

The  lines  of  constant  i  were  constrained  to  be  always  normal  to  the  wall.  This  is 
necessary  within  the  boundary  layer  for  stability  analyses.  The  distribution  of  points 
on  the  body  surface  was  specified,  and  the  distribution  on  the  outer  boundary  was 
determined  by  the  wall-normal  requirement.  For  blunt  nose  grids,  the  surface  spacing 
was  uniform  on  the  spherical  portion  of  the  nose,  and  a  form  of  exponential  stretching- 
using  two  segments  was  employed  to  expand  the  spacing  downstream. 

The  lines  of  constant  j  were  constructed  by  interpolating  between  spacings  spec¬ 
ified  at  the  upstream  and  downstream  boundaries.  Exponential  stretching  was  used 
to  cluster  points  within  the  boundary  layer.  An  experiment  was  conducted  in  which 
grid  points  were  clustered  at  both  the  wall  and  the  boundary-layer  edge.  A  mean 
flow  solution  was  analyzed  to  give  the  edge  location,  and  that  data  was  fed  along 
with  the  axial  point  distribution  into  the  commercial  software  GRIDGEN.  The  grid 
was  constructed  in  GRIDGEN  and  exported,  and  the  resulting  file  was  converted  to 
the  format  required  by  the  mean  flow  solver.  The  process  was  extremely  cumber¬ 
some,  and  no  significant  difference  was  seen  in  the  stability  results.  The  large  total 
number  of  grid  points  used  reduced  the  need  for  more  advanced  clustering  schemes. 

Figure  2.1  shows  examples  of  typical  grids  for  both  blunt  and  sharp  cones.  For 
clarity,  only  every  tenth  point  is  shown,  and  the  range  is  limited  to  the  nosetip  region. 

2.2.2  Mean  Flow  Solver 

A  modified  form  of  the  mean  flow  solver  DPLR2D  is  included  with  STABL. 
DPLR2D  solves  the  finite-volume  axisymmetric  Navier-Stokes  equations  for  real  gas 
flows  in  chemical  and  thermal  nonequilibrium  using  a  general  chemistry  model. 

DPLR2D  is  based  on  the  data-parallel  line  relaxation  method,  described  by 
Wright  et  al.  [25].  The  use  of  an  implicit  method  allows  for  rapid  convergence  of 
the  stiff  equation  set  that  results  from  the  vastly  different  length  scales  encountered 
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Figure  2.1:  Representative  grids  for  blunt  and  sharp  cone  cases.  The  range  is  limited 
to  the  nosetip  region,  and  only  every  tenth  grid  point  is  shown  for  clarity. 


in  high  Reynolds  number  flows.  DPLR2D  takes  advantage  of  the  structured  mesh 
to  provide  for  efficient  interprocessor  communication  without  explicit  domain  de¬ 
composition.  The  resulting  method  combines  the  parallelization  of  the  data-parallel 
lower-upper  relaxation  method  with  the  high  cell  aspect  ratio  convergence  prop¬ 
erties  of  the  Gauss-Seidel  line  relaxation  method.  DPLR2D  parallelizes  extremely 
efficiently,  allowing  for  very  large  problem  sizes  constrained  primarily  by  available 
memory. 

A  modified  form  of  Steger-Warming  flux  vector  splitting  is  used  to  give  an  accu¬ 
rate  solution  within  the  boundary  layer.  The  original  flux  vector  splitting  method 
proposed  by  Steger  and  Warming  [27]  was  developed  primarily  for  the  solution  of 
inviscid  flows.  MacCormack  and  Candler  [28]  showed  that  flux  splitting  causes  prob¬ 
lems  in  boundary  layers  through  either  abnormally  large  numerical  mixing  or  a  ficti¬ 
tious  pressure  gradient.  Their  modification  makes  the  method  much  more  accurate 
within  boundary  layers.  Second  order  accuracy  is  used  in  both  the  streamwise  and 
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wall-normal  directions,  but  a  limiter  switches  the  algorithm  to  first  order  accuracy 
near  the  shock. 

DPLR2D  has  excellent  convergence  properties.  The  highest  stable  Courant- 
Friedrichs-Lewy  (CFL)  number  begins  near  one  and  rapidly  increases  as  the  solution 
converges.  A  set  of  CFL  numbers  can  be  specified  in  the  input  file,  and  a  suggestion 
by  the  author  led  to  a  modification  to  allow  a  user  to  change  the  CFL  number  on- 
the-fly,  which  proved  helpful  when  working  on  a  time-sharing  system.  Stable  CFL 
numbers  of  at  least  10,000  are  attainable,  although  2,000  was  usually  the  maximum 
used  to  obtain  solutions  for  this  thesis.  The  maximum  stable  CFL  number  decreases 
with  increasing  mesh  density  and  aspect  ratio,  so  DPLR2D  offers  the  option  to  freeze 
a  portion  of  the  solution  from  the  stagnation  point  to  a  specified  distance  once  that 
portion  is  sufficiently  converged.  This  freezing  stops  the  iterations  in  the  specified 
area,  which  eliminates  the  highest  aspect-ratio  cells  from  the  solution  space,  allowing 
for  higher  CFL  numbers  and  faster  convergence  of  the  downstream  section.  Although 
this  option  is  helpful,  it  was  not  used  for  most  of  the  cases  presented  in  this  thesis 
since  it  was  easier  to  let  the  solution  run  overnight  unattended.  This  option  will 
probably  be  more  beneficial  for  design  purposes  when  the  amount  of  analysis  time 
needed  for  any  given  configuration  is  less. 

Boundary  conditions  are  specified  at  the  freestream,  outflow,  and  wall.  At  the 
wall,  a  no-slip  condition  is  imposed  along  with  requirements  for  zero  wall-normal 
pressure  gradients  and  mass  concentration  gradients.  Either  an  adiabatic  wall  or 
a  specified  wall  temperature  was  assumed,  as  indicated  for  each  individual  case. 
The  Mach  number,  static  density,  static  translational  and  vibrational  temperatures, 
and  species  mass  fractions  are  specified  in  the  freestream.  The  inflow  and  outflow 
boundaries  are  supersonic  except  in  the  boundary  layer.  A  zero  gradient  condition 
is  applied  at  the  outflow  boundary.  Zero  gradient  is  also  assumed  along  the  outer 
domain  boundary  to  ensure  relaxation  of  the  flow  variables. 

Real  gas  effects  are  modeled  following  the  method  of  Candler  and  MacCor- 
mack  [29].  These  effects  include  non-equilibrium  chemical  reactions  and  thermal 
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Figure  2.2:  Typical  RMS  residual  and  CFL  behavior  for  the  mean-flow  computations 
described  in  this  thesis. 


non-equilibrium  and  have  been  shown  to  be  critical  factors  in  obtaining  accurate 
predictions  of  boundary-layer  stability  and  transition  for  high  enthalpy  flows  [21]. 
For  the  cases  presented  here,  a  five-species  air  model  consisting  of  N2,  02,  N,  O, 
and  NO  was  used  with  the  standard  freestream  mass  fractions  of  76.7%  N2  and 
23.3%  02  (79%  N2  and  21%  02  by  volume)  [24],  Thermal  equilibrium  was  assumed 
in  the  freestream  except  where  otherwise  specified.  Section  3.2.1  gives  a  detailed 
discussion  of  the  transport  properties  used.  Stoke’s  hypothesis  is  used  for  the  bulk 
viscosity.  Equilibrium  constants  and  other  properties  are  taken  from  the  NASA  ther¬ 
modynamic  data  file  for  the  CEA  program  [30].  The  stability  solver  also  used  this 
chemistry  modeling. 

The  residuals  were  checked  in  two  stages  to  ensure  iterative  convergence.  DPLR 
prints  a  convergence  file  that  includes  the  RMS  residual  based  on  the  density.  That 
value  was  monitored  during  the  run,  and  a  typical  convergence  history  can  be  seen 
in  Figure  2.2.  A  residual  drop  of  approximately  14  orders  of  magnitude  is  present. 
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The  spikes  within  the  first  5000  iterations  are  thought  to  be  due  to  fluctuations  in 
the  location  of  the  shock  relative  to  the  grid  lines.  The  CFL  number  started  small 
but  was  rapidly  increased  by  five  orders  of  magnitude.  As  a  second  check  on  the 
convergence,  after  the  RMS  had  flattened,  contour  plots  of  the  residual  were  plotted 
using  an  exponential  scale  to  look  for  any  local  spots  of  higher  residual  and  to  ensure 
that  the  solution  was  fully  converged  on  the  aft  end  of  the  body. 


2.2.3  Stability  Solver 


The  code  PSE-Chem  is  used  to  analyze  the  stability  of  the  laminar  mean  flow 
profiles.  PSE-Chem  uses  the  linear  PSE  coupled  with  the  eN  method  to  produce 
transition  predictions.  Since  the  PSE  constitute  an  initial-boundary-value  problem, 
the  wavenumber  and  eigenfunction  at  the  starting  location  are  obtained  from  LST. 
A  summary  of  the  relevant  equations  and  their  derivation  will  be  given  here  for 
familiarization.  The  full  derivation  used  in  PSE-Chem  is  given  by  Johnson  [23]. 

The  stability  equations  are  obtained  from  the  Navier-Stokes  equations  by  first 
decomposing  the  instantaneous  flow  into  a  mean  and  a  fluctuating  component, 
q  =  q  +  q',  where  q  is  any  flow  variable.  This  decomposed  form  is  substituted  into 
the  Navier-Stokes  equations,  and  the  mean  flow  equation  is  subtracted,  resulting  in 
the  disturbance  equation  given  in  Equation  2.1,  where  <f>  is  the  vector  of  disturbance 
quantities  given  in  Equation  2.2. 


dcf)  d(j>  d(j>  d(j>  _  d2(f)  T j.  d2(f) 

Tm+Adi  +  %  +  cvFD,l,+  VlId^  +  V”W 


+  +  v. 

oz2 


d2<f> 
xy  dxdy 


+  Vm.*L  +  v„*i- 

oxo z  oyoz 


+  Fr 


<t>  =  {p\:P2:  ■■■:  Pus,  U' ,  V\  W' ,  T\  T'v) 


i\T 


0  (2.1) 
(2.2) 


The  terms  in  the  Jacobian  matrices  T,  A,  B, . . . ,  Vyz  depend  only  on  the  known  mean- 
flow  variables  and  their  derivatives,  and  all  of  the  terms  that  vary  nonlinearly  with 
(j)  are  contained  in  the  vector  Fn. 
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The  disturbance  quantities  are  assumed  to  be  traveling  waves  of  the  form  given 
by  Equation  2.3. 

(2.3) 


Here  £  and  r)  are  the  body-tangential  and  body-normal  computational  coordinates, 
k  is  the  real  spanwise  wavenumber,  z  is  the  distance  in  the  spanwise  direction,  and 
uj  is  the  real  frequency. 

When  Equation  2.3  is  substituted  into  Equation  2.1,  the  derivatives  are  cal¬ 
culated,  and  the  equation  is  converted  to  computational  coordinates,  the  result  is 
Equation  2.4,  where  the  Jacobians  are  given  in  Reference  [23]. 


vv 


P±  +  v  ^  i  r ' 

drf  ^  d£dr) 


(2.4) 


X  is  decomposed  using  Equations  2.5  and  2.6, 


X  =  ip(t,r])A(0  (2.5) 

A(£)  =  eie{i)  (2.6) 

where  d6/d £  =  «(£),  a  is  the  body-parallel  wavenumber  in  computational  coordi¬ 
nates,  and  ip  is  the  shape  function  vector.  This  allows  for  all  of  the  ellipticity  in  the 
wave  function  to  be  retained  while  parabolizing  only  the  shape  factor.  After  further 
substitutions  and  rearranging,  the  result  is  the  disturbance  equations  in  computa¬ 
tional  coordinates, 


A  /  7^  -  d2ip  -  d2ip  -  d2ip 

Di,  +  A—  +  B  —  +  Vt(—  +  Vm—  +  Va— 


+  Fn  =  0 


(2.7) 


To  arrive  at  the  PSE,  the  terms  and  V^n  are  neglected.  The  distur¬ 

bances  are  assumed  to  be  small,  allowing  the  nonlinear  terms  that  comprise  Fn  to 
be  neglected  as  well.  The  resulting  terms  give  the  linear  PSE, 


—  ~dib  ''dip 
Dip  +  A— +  -B— —  +  Vr 
a£  dr) 


d2ip 

vvW 


(2.8) 


The  PSE  constitute  an  initial  boundary  value  problem,  and  the  marching  pro¬ 
cedure  requires  an  initial  wavenumber  solution,  which  is  obtained  from  the  linear 
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stability  equations.  These  are  obtained  by  assuming  a  quasi-parallel  flow,  which 
means  that  ip  =  ip(rj)  and  da/dt;  =  0.  Equation  2.9  gives  the  resulting  linear  stabil¬ 
ity  equations. 


-  ^dip  -  d2ip 
DtP  +  B1f  +  Vvri^r  =  0 


(2.9) 


dp  '  m  dr/2 

These  equations  are  solved  using  the  global  and  local  two-step  procedure  of  Ma¬ 
lik  [31], 

In  the  global  procedure,  the  approximation  that  a2  =  0  is  made.  This  results 
in  a  generalized  eigenvalue  problem  for  the  complex  wavenumber  a.  Second-order 
central  differences  are  used  to  discretize  the  T]  derivatives.  The  boundary  condi¬ 
tions  of  zero  velocity  and  temperature  disturbances  and  zero  disturbance  pressure 
gradient  are  imposed  at  the  wall.  Zero  disturbances  are  imposed  for  all  variables 
at  the  outer  domain.  The  generalized  eigenvalue  problem  is  solved  using  the  LZ 
algorithm  [32],  The  result  is  a  spectrum  of  approximate  eigenvalue  guesses  for  the 
complex  wavenumber  a. 

The  local  procedure  is  then  used  to  refine  the  wavenumber  guesses  produced 
by  the  global  method.  Since  the  global  process  is  computationally  expensive,  it 
is  run  for  a  small  number  of  frequencies  and  PSE-Chem  interpolates  those  results 
to  intermediate  frequencies.  Discretization  is  accomplished  through  fourth-order 
central  differences  at  the  interior  points  and  second-order  central  differences  at  the 
boundaries.  A  non-homogeneous  boundary  condition  must  be  imposed  to  prevent 
arriving  at  the  trivial  solution  for  Equation  2.9.  The  u  =  0  condition  is  replaced  with 
a  normalization  condition  for  the  density  fluctuation  p  —  pnorm  =  p^.  The  resulting 
system  is  a  block  pentadiagonal  matrix  system  containing  terms  in  both  a  and  a2. 
Using  the  initial  guess  provided  by  the  global  method,  iterations  are  performed  on  a 
using  the  secant  method  until  is  reduced  to  zero.  This  results  in  the  complex 

wavenumber  a  and  the  disturbance  eigenfunction  ip  at  the  starting  location. 

The  PSE  are  solved  using  the  method  of  Herbert  [13].  The  variables  a  and  ip  are 
updated  simultaneously  at  each  marching  step.  The  boundary  conditions  of  LST  are 
used,  but  without  the  normalization  of  p.  The  disturbance  kinetic  energy  is  used  to 
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normalize  the  integrals  used  in  the  wavenumber  updates.  Iterations  are  performed 
until  the  wavenumber  has  converged  to  within  some  tolerance,  and  that  solution  is 
used  as  the  starting  guess  for  the  next  marching  station.  Further  details  of  the  PSE 
solution  procedure  implemented  in  PSE-Chem  are  given  in  Reference  [23]. 


Transition  prediction  is  accomplished  using  the  eN  method.  The  value  eN ,  where 
N  is  given  by  Equation  2.10,  represents  the  total  growth  factor  of  a  small  amplitude 
initial  disturbance. 


N  =  j*  a  di 

(2.10) 

For  the  PSE,  the  disturbance  growth  rate  a  is  obtained  from  Equations  2.11  and  2.12. 

Here  the  domain  of  integration,  fi,  extends  in  the  body-normal  direction  from  the 

surface  to  the  outer  boundary. 

1  dE 

(2,11) 

a  =  —otn  H - — 

?  2  E  di 

E  =  f  p(\u\2  +  \v\2  +  \w\2)di 

(2.12) 

n 


This  differs  from  LST,  where  only  the  contribution  of  -cq  is  used  to  calculate  N 
factors,  as  opposed  to  including  the  change  in  the  disturbance  shape  factor.  The 
disturbance  kinetic  energy  is  only  one  option  available  for  calculating  the  amplitude. 
Johnson  [26]  showed  that  the  use  of  the  massflux  can  produce  slightly  different 
results,  but  the  ability  to  perform  the  calculations  either  way  is  useful  for  comparisons 
with  other  results. 

Mack  [12]  showed  that  the  disturbances  with  the  form  given  in  Equation  2.3  can 
have  several  modes.  The  first  mode  is  analogous  to  the  Tollmien-Schlichting  waves 
of  incompressible  flow.  It  is  damped  by  wall  cooling  and  is  most  amplified  when  it 
is  at  an  oblique  angle.  The  second  mode  can  be  thought  of  as  a  trapped  acoustical 
wave.  It  is  amplified  by  wall  cooling,  and  it  is  most  amplified  when  it  is  2D.  Higher 
modes  exist,  but  they  always  have  lower  amplification  rates  than  the  first  or  second 
mode.  Only  the  first  and  second  modes  will  be  considered  in  this  thesis. 

The  stability  solver  uses  a  different  grid  than  the  mean  flow  solver.  Before  be¬ 
ginning  the  global,  local,  and  PSE  marching  steps,  the  mean  flow  solution  at  each 
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streamwise  location  is  linearly  interpolated  onto  a  grid  with  parameters  specified  in 
the  PSE-Chem  input  file.  This  grid  can  have  standard  hyperbolic  or  exponential 
stretching,  or  it  can  cluster  points  at  both  the  wall  and  edge  using  either  hyper¬ 
bolic  or  cosine-exponential  stretching.  This  interpolation  will  slightly  decrease  the 
accuracy  of  the  mean-flow  variables,  but  the  effect  of  this  change  was  not  examined. 
No  significant  difference  was  found  in  the  results  obtained  with  the  various  available 
stretching  methods. 

PSE-Chem  is  designed  for  parallel  processing  using  the  MPI  message  passing 
libraries.  The  global  solution  procedure  divides  the  frequencies  to  be  searched  among 
the  available  processors.  The  local  procedure  divides  the  guesses  to  be  refined  for  a 
given  frequency  among  the  processors.  The  linear  PSE  marching  does  not  parallelize 
as  well,  but  this  step  is  normally  much  faster  than  the  global  or  local  steps. 

The  majority  of  the  results  presented  in  this  thesis  were  obtained  using  Perl 
scripts  to  automatically  run  many  cases  sequentially,  label  the  cases  logically,  and 
save  the  relevant  results.  Appendix  A  gives  an  example  of  a  script  to  calculate  N 
factors  for  combinations  of  starting  location,  co,  and  (3.  These  scripts  use  functions 
from  the  scripting  library  recently  added  to  STABL.  This  is  a  powerful  and  robust 
methodology,  but  it  is  inefficient.  The  computations  performed  for  this  thesis  gener¬ 
ally  used  a  wide,  but  fairly  coarse,  band  of  parameters.  This  ensured  that  the  most 
unstable  modes  were  captured,  but  the  resolution  near  the  most  unstable  modes  was 
not  always  as  good  as  desired.  To  increase  the  resolution  with  this  system  without 
resulting  in  an  unreasonable  number  of  cases  would  require  either  better  knowledge 
of  the  expected  results  to  reduce  the  parameter  bandwidth  or  the  development  of 
some  sort  of  adaptive  refinement  algorithm.  With  the  current  system,  the  best  solu¬ 
tion  would  probably  be  to  run  a  wide,  coarse  parameter  matrix,  analyze  the  results, 
and  then  run  a  second  set  of  cases  with  a  finer  parameter  matrix  centered  on  the 
critical  modes. 
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2.3  Equipment  Used 

The  computations  described  in  this  thesis  were  performed  on  one  of  two  systems. 
A  Dell  Precision  650  workstation  with  two  3.06  GHz  Intel  Xeon  processors  and  2 
GB  of  RAM  was  used  for  some  of  the  computations  and  all  of  the  post-processing. 
In  addition,  access  to  the  ‘Rogue’  Beowulf  cluster  was  provided  by  Sandia  National 
Laboratories.  This  cluster  has  247  nodes  with  two  3.06  GHz  Intel  Xeon  processors 
on  each  node,  connected  by  2  Gb/s  Myrinet  connections,  and  160  nodes  with  two  2.4 
GHz  Intel  Xeon  processors  on  each  node,  connected  by  100  Mb/s  Ethernet  connec¬ 
tions.  All  nodes  on  Rogue  had  2  GB  of  RAM.  Both  the  local  workstation  and  Rogue 
used  the  Red  Hat  9  Linux  operating  system.  The  workstation  had  a  2.4.20-31.9smp 
kernel,  and  Rogue  had  a  2.4.20-30.9smp  kernel  with  GPFS  patches  applied.  The 
codes  were  compiled  using  the  Portland  Group  pgf  90  and  pgcc  compilers.  Message 
passing  was  accomplished  using  the  MPICH  MPI  libraries. 

Both  mean  flow  and  stability  calculations  were  performed  on  both  systems.  No 
significant  difference  was  seen  between  the  results  on  either  system.  Solutions  with 
150,000  cells,  which  were  typical  of  the  size  used  in  this  thesis,  could  be  obtained 
in  approximately  four  hours  on  sixteen  processors,  and  the  solution  time  scaled 
approximately  linearly  with  the  number  of  processors.  The  residual  typically  dropped 
by  12-14  orders  of  magnitude.  If  necessary,  the  computational  time  required  could  be 
significantly  decreased  by  freezing  the  solution  and  increasing  the  CFL  number  more 
aggressively.  The  global  and  local  solutions  each  took  approximately  one  minute  per 
frequency  per  processor.  The  time  required  for  sweeps  over  multiple  combinations  of 
starting  location,  frequency,  and  spanwise  wavenumber,  as  shown  in  Chapter  4,  was 
approximately  16  hours  on  two  2.4  GHz  processors.  A  more  advanced  script  using 
an  adaptive  refinement  scheme  could  significantly  decrease  the  run  time  if  necessary. 

Appendix  B  provides  a  detailed  description  of  the  methods  of  operation  used  for 
STABL  as  well  as  the  lessons  learned  along  the  way. 


19 


3.  Second  Mode  Verification:  Stetson’s  Blunt  Cone 

Chapters  3  and  4  will  present  the  results  of  computations  that  were  performed  to 
verify  and  validate  the  ability  of  STABL  to  accurately  compute  instabilities  and 
predict  transition.  All  of  the  cases  are  for  wind-tunnel  conditions.  Wind-tunnel 
cases  were  chosen  for  verification  because  a  large  body  of  previous  computations 
exists  and  to  match  the  general  conditions  of  the  Rufer  wind-tunnel  experiments, 
which  have  not  been  previously  analyzed.  The  Rufer  results  will  be  presented  in 
Chapter  5.  Since  it  is  not  always  possible  to  predict  which  is  the  dominant  mode, 
computations  for  both  the  first  and  second  modes  were  verified. 

The  second  mode  calculations  in  STABL  were  verified  by  benchmarking  against 
other  codes  since  accurate  experimental  data  were  not  available  for  validation.  This 
exercise  uncovered  a  number  of  issues,  including  bugs,  alternative  algorithms,  and 
different  transport  property  models.  These  discoveries  led  to  changes  that  improved 
the  quality  of  STABL.  However,  this  was  not  an  ideal  situation,  as  many  of  the 
details  of  the  other  calculations  are  unknown.  The  methods  used  are  not  identical, 
as  STABL  uses  the  PSE  and  is  designed  for  chemically  reacting  flows,  whereas  all  of 
the  benchmark  codes  use  LST  for  perfect  gas  air  flows.  In  addition,  unknown  errors 
may  exist  in  the  various  codes  or  their  usage.  These  could  account  for  the  remaining 
differences  in  the  results. 

3.1  Experimental  Conditions 

The  conditions  chosen  for  the  exercise  were  those  of  the  blunt  cone  experiments 
of  Stetson  et  al.  [33]  with  the  3.81  mm  nose  radius.  Stetson  performed  hot-wire 
measurements  in  the  boundary  layer  of  a  sphere-cone  at  0°  angle  of  attack  and  Mach 
8  in  Tunnel  B  at  the  Arnold  Engineering  and  Development  Center  (AEDC).  The  test 
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Table  3.1:  Test  conditions  for  Stetson  et  al.  [33]  blunt  cone  experiment  for  all  com¬ 
putations  presented  here. 


Model  Specifications 

Half  angle 

7° 

Length  (m) 

1.016 

Base  diameter  (m) 

0.2495 

Nose  radius  (mm) 

3.81 

Wall  Temperature 

Adiabatic 

Freestream  Conditions 

Fluid 

Air 

Moo 

7.99 

Po  (Pa) 

4.0xl06 

To  (K) 

750 

Poo  (Pa) 

410 

Too  (K) 

54 

Poo  (kg/m3) 

0.0027 

Reoo/ft 

2.5xl06 

conditions  for  all  Stetson  computations  in  this  thesis  are  summarized  in  Table  3.1. 
The  freestream  static  temperature  and  density  are  required  inputs  for  STABL,  and 
they  were  calculated  from  the  specified  flow  conditions  using  the  isentropic  relations 
and  the  perfect  gas  law.  Further  discussion  of  the  input  conditions  will  be  given  in 
Section  3.2.2. 

The  Stetson  et  al.  experiment  has  been  analyzed  by  numerous  other  researchers 
using  many  different  codes,  which  makes  it  an  ideal  benchmark  for  wind-tunnel  cases. 
All  of  the  other  researchers  used  LST-based  codes  for  stability  calculations,  as  op¬ 
posed  to  the  PSE  used  in  STABL.  Malik  et  al.  [34]  performed  the  first  computations. 
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To  compute  the  mean  flow,  they  used  a  full  Navier-Stokes  solution  for  flow  around 
a  sphere  to  provide  an  initial  data  plane  for  a  Parabolized  Navier-Stokes  (PNS)  cal¬ 
culation  on  the  cone  frustum.  Esfahanian  [35]  performed  Thin-Layer  Navier-Stokes 
calculations  using  very  fine  grid  spacing  by  the  standards  of  the  day.  Kufner  et 
al.  [36]  performed  a  careful  analysis  of  the  sensitivity  of  the  stability  results  to  vari¬ 
ous  assumptions  made  in  computing  the  mean  flow.  Results  obtained  by  Stilla  using 
a  second-order  boundary-layer  solver  are  included  in  Kufner  et  al.  [36],  as  well  as 
in  Reference  [37].  Rosenboom  et  al.  [38]  used  Kufner’s  results  as  a  basis  for  exam¬ 
ining  the  effect  of  nose-tip  bluntness  on  stability.  More  recently,  Schneider  [4,  39] 
performed  a  detailed  reassessment  of  the  Stetson  experimental  data  and  performed 
stability  analyses  using  the  mean  flow  generated  by  a  boundary-layer  solver.  Es¬ 
fahanian  [40]  performed  additional  computations  using  Stetson’s  conditions  to  test 
the  accuracy  of  PNS  mean  flow  solutions  for  stability  applications.  Lyttle  et  al.  [41] 
used  a  DPLR  solver  to  analyze  Stetson’s  case  with  a  particular  emphasis  on  the  wall 
temperature  distribution.  Zhong  [42]  performed  LST  calculations  as  part  of  an  effort 
to  analyze  the  flow  using  Direct  Numerical  Simulation. 

Schneider  [39]  has  shown  that  there  are  a  number  of  open  issues  associated  with 
this  experiment.  The  axial  station  s/rn  =  175  is  used  for  the  vast  majority  of  the 
computations  presented  in  the  literature.  However,  Stetson’s  data  show  considerable 
nonlinearity  at  this  station  which  will  not  be  captured  by  a  LST  analysis.  Overshoots 
are  present  in  the  pitot-pressure  profiles,  which  are  thought  to  be  due  to  interference 
from  the  probe  and  have  not  appeared  in  any  computation.  Most  researchers  have 
assumed  the  wall  is  adiabatic,  but  the  measured  surface  temperature  is  about  20% 
below  the  adiabatic  temperature.  The  present  work  is  not  an  attempt  to  resolve  the 
discrepancies  between  computations  and  experiments.  Rather,  the  intent  is  to  repli¬ 
cate  the  computational  results  obtained  with  other  codes  to  partially  verify  STABL. 
The  large  body  of  independent  computational  data  makes  this  case  ideal  for  verifi¬ 
cation,  but  these  unresolved  issues  make  complete  validation  from  the  experimental 
data  impossible. 
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Figure  3.1:  Amplification  rate  for  Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  =  3.81 
mm  as  computed  by  other  researchers.  Symbols  are  for  identification  only. 


Most  of  the  researchers  mentioned  above  published  plots  of  the  amplification 
rate  at  s/rn  =  175.  Figure  3.1  shows  these  results  on  a  single  plot.  The  data  was 
obtained  by  digitizing  the  figures  in  the  references  using  the  Un-Scan-It  software 
package.  The  figures  were  electronically  scanned  to  high  resolution,  and  then  the 
Un-Scan-It  software  was  used  to  convert  the  graphical  data  to  numerical  data  using 
an  algorithm  that  interpolates  based  on  the  pixel  coordinates  of  the  line  and  the 
axes.  The  accuracy  of  the  digitization  is  estimated  to  be  on  the  order  of  the  line 
width. 

The  computations  in  Figure  3.1  agree  much  better  with  each  other  than  they 
do  with  the  experimental  data.  This  is  most  likely  due  to  non-linearity  and  the 
other  issues  discussed  above.  The  average  peak  amplification  rate  based  on  the  six 
computations  plotted  is  14.4/m  at  a  frequency  of  129  kHz. 
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3.2  Issues  Addressed 

A  number  of  issues  were  addressed  to  try  to  improve  the  agreement  between  the 
stability  results  of  STABL  and  those  of  other  researchers.  Results  incorporating  the 
effects  of  the  lessons  learned  from  all  of  these  issues  are  presented  in  Section  3.3.  The 
following  sections  are  intended  to  illustrate  the  incremental  effect  of  each  individual 
issue.  Therefore,  each  section  should  be  viewed  as  distinct  from  the  others,  rather 
than  as  a  sequential  progression  of  increasingly  accurate  solutions. 

3.2.1  Laminar  Transport  Properties 

The  laminar  transport  properties  are  an  essential  component  of  any  numerical 
model.  These  properties  are  the  viscosity  for  momentum  transport,  the  thermal 
conductivity  for  energy  transport,  and  the  diffusion  coefficient  for  mass  transport. 
Accepted  models  exist  for  usage  within  normal  ranges  of  temperature  and  pressure. 
However,  there  is  considerable  subjectivity  in  determining  the  best  model  for  very 
high  or  low  temperatures,  and  the  best  choice  for  one  extreme  is  generally  not  the 
best  choice  for  the  other.  Lyttle  and  Reed  [43]  showed  for  high  temperatures  that 
the  transport  properties  employed  can  have  a  large  impact  on  hypersonic  stability 
results. 

The  components  of  STABL  were  designed  primarily  for  high  enthalpy  air  flows, 
such  as  those  found  in  hypersonic  flight  or  in  shock  tunnels.  These  flows  are  charac¬ 
terized  by  extremely  high  temperatures,  causing  dissociation,  ionization,  and  chem¬ 
ical  and  thermal  non-equilibrium,  phenomena  known  collectively  as  real-gas  effects. 
Much  more  complicated  flow  models  are  required  to  accurately  model  these  real  gas 
effects  than  for  situations  where  perfect  gas  behavior  can  be  assumed. 

Conventional  hypersonic  wind  tunnels  obtain  high  Mach  numbers  in  part  by 
lowering  the  static  temperature  in  the  test  section  to  just  above  the  condensation 
point  of  nitrogen.  This  lowers  the  speed  of  sound,  creating  hypersonic  flow  at  lower 
freestream  velocities  than  would  be  found  in  atmospheric  flight  at  the  same  Mach 
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number  [1],  Tunnel  B  at  AEDC  and  the  Purdue  Mach-6  Ludwieg  Tube  (M6LT)  both 
have  freestream  static  temperatures  below  60  K  for  standard  operating  conditions. 
For  slender  bodies  with  highly  oblique  shocks,  boundary-layer  edge  temperatures  of 
60-100  K  are  typical.  For  the  STABL  suite  to  be  valid  in  the  low  enthalpy  regime, 
transport  property  models  appropriate  to  that  regime  must  be  employed. 

All  of  the  computations  cited  previously  were  performed  with  a  Sutherland  vis¬ 
cosity  law.  Not  all  of  the  references  provide  details  about  the  thermal  conductivity 
modeling,  but  those  that  did  specify  the  method  assumed  a  perfect  gas  with  a  con¬ 
stant  Prandtl  number.  It  is  probable  that  those  who  did  not  specify  their  method 
used  this  technique  as  well.  Because  all  of  the  computations  in  the  literature  treated 
the  air  as  a  single  gas  rather  than  a  mixture,  none  included  the  effects  of  diffusion. 
Due  to  the  low  temperatures  involved,  the  production  of  monatomic  species  is  neg¬ 
ligible.  The  binary  diffusion  model  should  therefore  not  be  a  significant  factor,  and 
its  accuracy  was  not  investigated. 


Viscosity 

STABL  originally  used  Blottner’s  [44]  viscosity  model  for  reacting  flows  to  de¬ 
termine  individual  species  viscosities  [23].  This  model  uses  empirical  curve  fits  of 
the  form  given  by  Equation  3.1  to  compute  individual  viscosities  in  kg/m-s  for  each 
species  in  the  mixture. 

fjla  =  0.1  e(AJn(T)+Bs)ln(T)+Cs  (3.1) 

The  viscosity  of  the  mixture  is  found  by  Wilke’s  semi-empirical  mixing  law,  given  in 
Equations  3. 2-3. 4. 
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Table  3.2:  Blottner’s  Viscosity  Curve  Fit  Coefficients 


Coefficient 

High  Temp 

Low  Temp 

A  n2 

0.0268142 

0.0089993 

Bv2 

0.3177838 

0.6039338 

C  n2 

-11.3155513 

-12.4453814 

Ao2 

0.044929 

-0.0255541 

Bo2 

-0.0826158 

1.0503525 

Co2 

-9.2019475 

-13.7080219 

r=  1 


(1+< 


(3,4) 


Two  sets  of  the  coefficients  As,  Bs,  and  Cs  are  given  in  Table  3.2.  The  high 
temperature  set  was  published  in  Blottner  et  al.  [44]  and  is  valid  over  the  tempera¬ 
ture  range  1000  K-30000  K.  The  low  temperature  set  represents  unpublished  data 
provided  by  Blottner  (personal  communication,  October  2004)  and  is  valid  over  the 
range  100  K-10000  K.  Values  for  the  species  NO,  N,  and  O  are  given  in  Reference  [44] 
and  are  the  same  for  both  temperature  ranges.  These  species  have  mass  fractions 
many  orders  of  magnitude  below  those  of  N2  and  02  at  wind  tunnel  temperatures, 
and  their  viscosity  modeling  was  not  investigated. 

The  Sutherland  law  is  the  method  most  commonly  used  in  CFD  codes.  This 
model  is  based  on  the  kinetic  theory  of  gases,  and  it  is  known  to  be  a  good  ap¬ 
proximation  in  the  temperature  range  100  K-2000  K  [15].  Its  form  is  given  by 
Equation  3.5,  where  S  is  an  effective  temperature  called  the  Sutherland  constant 
and  /iref  and  Tref  are  reference  values  [45]. 


F  _  (  T  \  2  Tref  +  S 
Hre /  \  Tre f  )  T  +  S 


(3.5) 
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Figure  3.2:  Viscosity  of  air  computed  using  various  methods  compared  with  experi¬ 
mental  data. 


Figure  3.2  shows  a  comparison  of  the  two  methods  given  in  Table  3.2  and  the 
Sutherland  law  plotted  with  the  experimental  data  of  Grieser  and  Goldthwaite  [46] 
and  Matthews  et  al.  [47].  Figure  3.3  shows  a  detail  of  the  same  plot  at  the  low 
temperature  limit.  The  viscosity  reported  by  Grieser  and  Goldthwaite  [46]  near  50 
K  is  two  orders  of  magnitude  smaller  than  the  next  closest  value.  This  was  based  on 
one  experimental  measurement,  and  it  illustrates  the  uncertainty  that  exists  at  very 
low  temperatures.  The  Blottner  high  temperature  model  is  seen  to  differ  greatly  from 
the  experimental  data  at  low  temperatures.  At  80  K,  which  is  approximately  the  edge 
temperature  for  the  Stetson  case,  the  Blottner  model  overpredicts  the  viscosity  by 
70%.  The  Sutherland  law  diverges  from  the  empirical  models  above  approximately 
1900  K,  as  is  expected  based  on  the  approximations  made  in  its  derivation  [45]. 

Based  on  this  analysis,  STABL  was  modified  to  incorporate  different  viscosity 
models  for  each  temperature  regime.  STABL  now  uses  the  Sutherland  law  for  low 
temperatures,  the  high-temperature  Blottner  law  for  high  temperatures,  and  the  low- 
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Figure  3.3:  Viscosity  models  and  experimental  data  for  air  at  low  temperatures. 


temperature  Blottner  law  for  the  intermediate  range.  Blending  functions  are  used  to 
ensure  a  smooth  transition  between  the  laws.  The  resulting  method  is  shown  as  the 
STABL  Blended  curve  in  Figures  3.2  and  3.3.  Agreement  with  the  experimental  data 
and  the  form  of  the  Sutherland  Law  used  in  the  other  computations  is  considerably 
better. 

It  is  important  to  note  that  there  is  considerably  higher  uncertainty  in  the  exper¬ 
imental  data  below  200  K  [48].  Several  models  exist  for  this  low  temperature  regime. 
The  Chapman-Enskog  model  [49]  is  based  on  the  kinetic  theory  of  gases  and  em¬ 
ploys  collision  integrals  to  model  the  intermolecular  forces.  Keyes  [50]  and  Maitland 
and  Smith  [51]  both  proposed  empirical  curve  fits,  and  Mack  [12]  proposed  a  model 
hereafter  referred  to  as  the  “Linear  Sutherland  Law”  that  is  a  linear  extrapolation 
of  the  Sutherland  Law  below  the  Sutherland  constant.  Because  the  experimental 
data  is  so  limited  for  this  temperature  range,  there  is  no  clear-cut  choice  for  the  best 
model.  The  Sutherland  law  was  used  for  all  computations  in  this  work  to  provide 
consistency  with  the  other  computations  referenced. 
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Thermal  Conductivity 


STABL  calculates  thermal  conductivities  for  both  translational  and  vibrational 
modes  [23].  It  uses  Eucken’s  relation,  given  in  Equations  3. 6-3. 8,  to  calculate  the 
translational  thermal  conductivity  of  each  species. 


—  /A  ^  2  ^vJrs  ^ 'v,rots  ^ 

3  R 


Cv 


rots 


Cv’trs  ~  2  Ms 

0  ,  monatomics 

jj-  ,  diatomics 
The  vibrational  thermal  conductivity  is  given  by  Equation  3.9. 

tie. 


k-vibs  —  l^s^v,viba  —  /A 


- vs 


dTv 


(3.6) 

(3.7) 

(3.8) 


(3.9) 


The  mixture  thermal  conductivities  are  computed  with  Wilke’s  semi-empirical  mix¬ 
ing  law,  where  ks  is  substituted  for  fis  in  Equation  3.2  and  Equations  3.3  and  3.4  are 
unchanged. 

The  other  researchers  that  have  computed  Stetson’s  blunt  cone  case  have  used 
perfect  gas  codes.  Those  that  did  specify  their  thermal  conductivity  method  used 
a  constant  Prandtl  number  of  0.72,  calculated  the  viscosity  using  Sutherland’s  law, 
assumed  a  constant  Cp.  and  calculated  the  thermal  conductivity  from  Equation  3.10. 

Pr  =  ^L  (3.10) 

k 


Figure  3.4  shows  the  variation  with  temperature  of  Prandtl  number  and  specific 
heat  at  constant  pressure,  for  air.  Data  of  Keyes  [52]  is  plotted  along  with  the  values 
used  by  STABL.  STABL  uses  a  Cp  that  is  constant  with  respect  to  temperature 
for  a  given  gas  mixture.  The  value  shown  was  that  calculated  by  the  Flow  Analysis 
module  of  the  GUI  for  a  mixture  of  0.767  N2  and  0.233  02.  The  Prandtl  number  was 
calculated  using  Equation  3.10  with  the  value  of  Cp  shown  in  Figure  3.4  and  viscosity 
and  thermal  conductivity  calculated  using  the  blended  model.  The  experimental 
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Figure  3.4:  Temperature  variation  of  Prandtl  number  and  specific  heat  at  constant 
pressure  for  air.  Data  is  taken  from  Keyes  [52], 


Prandtl  number  varies  from  0.77  at  100  K  to  0.68  at  500  K,  a  range  of  12%  within  the 
temperature  span  of  the  Stetson  flowfield.  The  specific  heat  is  reasonably  constant  at 
low  temperatures,  but  it  too  begins  to  vary  considerably,  differing  from  the  assumed 
constant  value  by  8%  at  750  K.  This  shows  that  the  common  assumptions  of  constant 
Cp  and  Pr  have  relatively  large  errors,  even  for  this  temperature  range.  The  impact 
of  these  assumptions  on  the  stability  results  should  be  examined. 

Figure  3.5  shows  a  comparison  of  the  thermal  conductivity  as  computed  by 
STABL  and  as  computed  using  a  constant  Prandtl  number  of  0.72.  The  curve 
marked  “STABL  Original”  uses  the  Blottner  high  temperature  viscosity  model,  the 
curve  marked  “STABL  Blended”  uses  the  blended  viscosity  model  described  earlier, 
and  the  constant  Prandtl  number  curve  uses  a  Sutherland  Law.  Experimental  data 
of  Taylor  and  Johnston  [53],  Keyes  [52],  and  Vines  [54]  are  also  presented.  Figure  3.6 
shows  a  detail  of  the  same  plot  for  the  low  temperature  regime.  The  STABL  blended 
method  agrees  very  well  at  low  temperatures  with  both  the  experimental  data  and 
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Figure  3.5:  Thermal  conductivity  for  air  computed  using  various  methods  compared 
with  experimental  data. 
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Figure  3.6:  Thermal  conductivity  models  and  experimental  data  at  low  tempera¬ 


tures. 


31 


Figure  3.7:  Velocity  profile  for  Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  =  3.81 
mm  with  original  viscosity  model  and  new  blended  model. 


the  constant  Prandtl  number  method.  The  use  of  the  blended  viscosity  law  causes 
the  error  in  the  calculated  thermal  conductivity  with  respect  to  the  experimental 
data  to  drop  from  70%  to  4%  at  80  K.  However,  all  three  of  the  methods  differ  sig¬ 
nificantly  from  the  experimental  data  at  higher  temperatures.  This  may  be  due  to 
the  importance  of  the  vibrational  energy  mode  at  higher  temperatures,  but  this  was 
not  investigated  further.  Given  the  close  agreement  at  high  temperatures  between 
the  constant  Prandtl  number  method  and  that  used  by  STABL,  this  difference  with 
the  experimental  data  should  not  be  a  source  of  difference  between  STABL  and  the 
computational  results  of  other  researchers.  However,  it  could  have  an  impact  on  the 
accuracy  of  actual  stability  predictions. 

Effects 

Figure  3.7  shows  the  effect  of  the  change  in  the  viscosity  model  on  the  boundary- 
layer  profile  at  s/rn  =  175.  The  blended  viscosity  model  makes  the  velocity  profile 
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Figure  3.8:  Second  derivative  of  the  velocity  in  the  normal  direction  for  Stetson’s 
blunt  cone  at  s/rn  =  175  with  rn  =  3.81  mm  with  the  original  and  blended  viscosity 
models. 


less  full,  and  the  agreement  with  the  data  of  Kufner  et  al.  [36]  improves  considerably. 
Figure  3.8  shows  the  effect  of  the  viscosity  model  on  the  second  derivative  of  the 
tangential  velocity  in  the  body-normal  direction.  The  shape  of  the  profile  agrees 
much  better  with  the  data  of  Kufner  et  al.  and  Esfahanian  [35].  The  blended  viscosity 
law  causes  the  magnitude  of  the  peak  of  the  second  derivative  curve  to  increase  by 
23%,  and  it  becomes  larger  than  that  of  Kufner  or  Esfahanian.  In  addition,  although 
agreement  is  improved  near  the  boundary-layer  edge,  there  are  still  small  differences 
in  the  profiles  in  this  region  critical  for  stability. 

Figure  3.9  shows  the  effect  of  the  viscosity  model  on  the  amplification  rate.  The 
amplification  curve  peak  is  shifted  from  (a;  kHz,— cq  l/m)  =  (136,13.4)  to  (132,17.3), 
which  is  an  increase  in  the  peak  amplification  rate  of  29%  and  a  decrease  of  3%  in 
the  frequency  at  which  that  occurs.  The  decrease  in  the  frequency  causes  better 
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Figure  3.9:  Amplification  rate  with  the  original  and  blended  viscosity  models  for 
Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  =  3.81  mm. 


agreement  with  the  other  researchers,  but  the  increase  in  the  amplification  rate 
causes  worse  agreement. 

3.2.2  Input  Conditions 

In  Reference  [33],  Stetson  et  al.  specify  the  quantities  M^,  P0 ,  T0,  P^,  and  unit 
Reynolds  number.  Discussions  with  Wayne  Hawkins  of  AEDC  (personal  commu¬ 
nication,  July  2005)  illustrated  how  these  numbers  were  obtained.  The  stagnation 
pressure  and  temperature  were  measured  directly  in  the  stilling  chamber.  The  Mach 
number  was  inferred  from  pitot  measurements.  The  freestream  static  quantities  were 
computed  using  a  Beattie-Bridgeman  correction  to  acount  for  real  gas  effects,  which 
was  less  than  a  1%  correction  at  these  conditions.  The  unit  Reynolds  number  is 
a  derived  quantity  that  depends  on  the  viscosity  model,  which  can  be  a  source  of 
uncertainty  at  the  freestream  temperature  of  54.4  K.  Table  3.3  shows  a  comparison 
of  the  unit  Reynolds  number  calculated  using  different  viscosity  laws  in  common  use 
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Table  3.3:  Unit  Reynolds  number  for  the  Stetson  et  al.  blunt  cone  experiment  cal¬ 
culated  using  various  viscosity  models. 


Viscosity  Model 

n  (  k9  ) 

Reoo/ft 

Sutherland 

3.546  x  10“6 

2.69  x  106 

Linear  Sutherland 

3.771  x  10“6 

2.53  x  106 

Keyes 

3.893  x  10“6 

2.45  x  106 

Chapman- Enskog 

3.788  x  10“6 

2.53  x  106 

for  the  low  temperature  regime.  Stetson  (personal  communication,  July  2005)  stated 
that  he  obtained  the  freestream  conditions  from  AEDC.  Hawkins  confirmed  that  the 
unit  Reynolds  number  was  calculated  using  a  linear  Sutherland  law. 

If  the  unit  Reynolds  number  is  used  to  determine  the  freestream  conditions,  the 
conditions  will  depend  on  the  viscosity  model  used.  All  of  the  other  researchers 
employed  a  standard  Sutherland  law  for  their  computations,  and  the  unit  Reynolds 
number  is  an  input  condition  for  many  CFD  codes.  To  truly  match  the  freestream 
conditions  of  Stetson’s  experiment  using  a  Sutherland  law  code,  a  unit  Reynolds 
number  of  2.69  x  106/ft  would  need  to  be  used,  rather  than  the  2.50  x  106/ft  specified 
in  Reference  [33].  It  seems  likely  that  at  least  some  of  the  other  computations  in  the 
literature  may  have  been  performed  using  incorrect  freestream  conditions. 

The  STABL  mean  flow  solver  uses  the  freestream  static  temperature,  freestream 
Mach  number,  and  freestream  static  density  as  input  conditions.  The  freestream 
static  density  is  0.02652^-  when  computed  from  the  stagnation  temperature  and 
pressure  using  the  perfect  gas  law  and  the  isentropic  relations.  If  a  Sutherland 
viscosity  law  is  used  with  a  unit  Reynolds  number  of  2.50  x  106/ft,  the  freestream 
density  becomes  0.02464^.  The  first  value  appears  to  be  more  physically  accurate, 
and  should  be  used  for  comparison  with  the  experimental  data.  However,  the  second 
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Figure  3.10:  Amplification  rates  computed  using  two  different  values  of  freestream 
density  for  Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  =  3.81  mm. 


number  is  most  likely  what  was  actually  used  in  the  other  researchers’  computations, 
and  thus  it  should  be  used  for  code  verification  with  the  current  benchmarks. 

Figure  3.10  shows  the  effect  of  changing  the  density  input.  Use  of  the  value  based 
on  the  unit  Reynolds  number  shifts  the  amplification  rate  curve  closer  to  the  other 
researchers’  results,  moving  the  peak  from  (132,17.9)  to  (127,16.6).  Figure  3.11  com¬ 
pares  N  factors  for  frequencies  near  100  kHz  calculated  with  each  of  the  freestream 
density  values.  Despite  the  difference  in  the  amplification  rates,  almost  no  difference 
is  evident  in  the  N  factors. 

A  similar  situation  was  described  by  Arnal  et  al.  [55].  When  they  attempted 
to  compute  stability  results  independently  by  two  groups,  the  use  of  two  different 
viscosity  models  between  the  codes  resulted  in  slightly  different  freestream  condi¬ 
tions.  Like  the  present  results,  this  caused  a  small,  but  not  neglible,  difference  in 
the  amplification  rate. 
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Figure  3.11:  N  factors  computed  using  two  different  values  of  freestream  density  for 
Stetson’s  blunt  cone  with  rn  =  3.81  mm. 


3.2.3  Boundary-Layer  Edge  Detection 

Many  aspects  of  stability  analysis  depend  on  an  accurate  determination  of  the 
location  of  the  boundary-layer  edge.  The  most  unstable  first  and  second  mode 
disturbances  occur  near  the  edge.  For  this  reason,  an  accurate  mean  flow  solution 
near  the  edge  is  vital,  and  PSE-Chem  contains  interpolation  algorithms  to  structure 
the  stability  grid  so  that  points  are  clustered  at  both  the  wall  and  edge.  In  addition, 
data  presented  in  the  literature  is  often  presented  in  terms  of  edge  quantities  such  as 
the  Reynolds  number  R,  but  the  edge  velocity,  temperature,  and  location  are  seldom 
specified. 

Unfortunately,  determining  the  location  of  the  boundary-layer  edge  is  not  triv¬ 
ial.  When  the  Prandtl  number  of  the  flow  is  not  unity,  the  viscous  and  thermal 
boundary-layers  will  have  different  thicknesses.  Valid  arguments  exist  for  using  each 
to  determine  the  “boundary-layer”  thickness.  In  addition,  a  blunt  body  introduces 
an  entropy  layer  that  must  be  considered.  Even  if  all  of  these  factors  are  considered, 
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a  parameter  must  be  specified  at  which  the  local  value  is  considered  to  be  the  edge. 
Choices  of  95%,  99%,  and  99.5%  are  all  common  in  the  literature. 

For  these  reasons,  STABL  avoids  the  use  of  the  edge  location  wherever  possible. 
Since  the  full  Navier-Stokes  equations  are  solved  over  the  entire  shock  layer,  an  edge 
parameter  is  not  used  in  determining  the  mean  flow.  The  edge  location  does  matter 
in  determining  the  grid  for  the  stability  analysis  when  edge  clustering  is  employed. 
However,  since  this  is  an  interpolation  from  the  existing  mean  flow  solution,  the  effect 
of  small  differences  in  the  edge  location  will  be  minimal.  The  wave  phase  speed  cr, 
given  by  Equation  3.11,  is  used  to  classify  instabilities  based  on  their  velocity  relative 
to  the  edge.  This  requires  an  accurate  determination  of  the  edge  location. 


STABL  uses  the  total  enthalpy  to  determine  the  boundary-layer  edge.  For  a 
Prandtl  number  less  than  unity,  the  thermal  boundary  layer  will  be  thicker  than 
the  viscous  boundary  layer.  In  addition,  the  total  enthalpy  is  constant  outside  of 
the  thermal  boundary  layer,  in  contrast  to  velocity,  which  may  vary  outside  of  the 
viscous  boundary  layer  due  to  pressure  gradient  effects. 

It  is  possible  for  a  flow  to  have  an  overshoot  in  its  total  enthalpy  profile,  as 
shown  in  Figure  3.12.  This  occurrence  leads  to  three  logical  possibilities  to  declare 
as  the  boundary-layer  edge,  as  indicated  by  the  horizontal  lines  in  Figure  3.12.  The 
lowest  line  corresponds  to  the  point  where  the  total  enthalpy  first  approaches  the 
freestream  value.  The  highest  line  is  where  the  total  enthalpy  begins  to  stabilize 
at  the  freestream  value.  The  middle  line  represents  the  peak  of  the  overshoot,  and 
could  be  seen  as  a  compromise  between  the  other  options. 

PSE-Chem  originally  used  the  lowest  line  as  the  boundary-layer  edge.  The  outer 
line  seems  to  be  a  much  better  choice,  as  shown  by  comparison  with  the  mass  flux 
profile  plotted  as  the  dashed  line  on  Figure  3.12.  Consequently,  the  PSE-Chem  algo¬ 
rithm  was  modified  to  choose  the  outer  line  as  the  edge  in  the  case  of  an  overshoot; 
it  still  handle  profiles  without  an  overshoot  correctly.  A  listing  of  the  algorithm 
proposed  by  the  author  is  given  in  Appendix  C. 
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Figure  3.12:  Representative  total  enthalpy  and  mass  flux  profiles. 


3.2.4  Grid  Generation 

The  requirement  for  accurate  second  derivatives  of  the  velocity  and  temperature 
profiles  throughout  the  boundary-layer  means  that  much  finer  grids  must  be  used 
for  stability  analysis  than  for  most  other  CFD  applications.  Although  results  could 
sometimes  be  obtained  with  coarser  grids,  80  points  in  the  boundary  layer  was  found 
to  be  an  approximate  minimum  to  obtain  grid  independent  results.  All  of  the  results 
presented  in  this  thesis  were  obtained  with  grids  containing  at  least  300  axial  points 
and  at  least  250  normal  points. 

Figure  3.13  shows  the  amplification  curve  computed  for  Stetson’s  blunt  cone  at 
s/rn  =  175  with  rn  =  3.81  mm  using  several  grids  with  varying  resolution.  The  five 
curves  shown  used  grids  with  106,  121,  140,  239,  and  276  points  in  the  boundary 
layer.  The  points  were  clustered  at  the  wall  and  exponentially  stretched  to  the  outer 
domain.  The  grids  marked  (b)  in  Figure  3.13  used  an  outer  grid  boundary  that  was 
very  carefully  constructed  to  match  the  shock  shape  over  the  entire  length  of  the 
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Figure  3.13:  Amplification  rate  for  Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  =  3.81 
mm  computed  using  mean  flows  with  varying  grid  resolutions. 


Figure  3.14:  N  factors  computed  using  mean  flows  with  varying  grid  resolutions  for 
Stetson’s  blunt  cone  with  rn  =  3.81  mm. 
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cone.  The  grids  marked  (a)  used  a  grid  that  matched  the  shape  only  reasonably  well 
and  is  representative  of  the  quality  of  grids  used  in  most  of  the  other  computations 
shown  in  this  thesis.  Small  differences  between  the  grids  are  present,  but  no  apparent 
trends  are  visible. 

Figure  3.14  shows  the  N  factors  computed  using  each  of  the  grids.  Again,  the 
differences  are  relatively  small.  The  grids  marked  (b)  produce  slightly  smaller  N 
factors  than  the  grids  marked  (a).  There  is  no  discernable  trend  with  the  total  grid 
size.  This  suggests  that  the  differences  are  due  to  the  sensitivity  of  the  stability 
calculations  to  small  grid  variations,  and  it  cannot  be  conclusively  said  that  one 
grid  shape  is  better  than  the  other.  The  scatter  of  approximately  5%  at  s  =  1.0 
m  may  represent  an  uncertainty  in  the  results  that  cannot  be  eliminated  through 
increasing  mesh  density.  In  any  case,  the  differences  are  slight  compared  with  the 
overall  uncertainty  associated  with  the  use  of  the  eN  method,  and  sufficient  mesh 
spacing  is  used  for  the  present  results. 

3.2.5  Wall  Temperature 

Most  previous  computational  comparisons  to  the  Stetson  experiment  have  as¬ 
sumed  an  adiabatic  wall,  but  Schneider  [4]  showed  that  there  was  actually  signifi¬ 
cant  heat  loss  through  the  sting.  Lyttle  et  al.  [41]  performed  stability  computations 
using  both  an  adiabatic  wall  and  wall  temperature  distribution  similar  to  the  ex¬ 
periment.  They  found  a  difference  of  approximately  15%  in  the  peak  amplification 
rate  at  s/rn  =  175  m,  but  the  difference  did  not  explain  the  differences  with  the 
experimental  amplification  rates.  All  of  the  Stetson  computations  in  this  thesis  were 
performed  assuming  an  adiabatic  wall  for  comparision  with  the  other  computations. 

Lyttle  et  al.  plot  their  calculated  adiabatic  wall  temperature  distribution  in  Fig¬ 
ure  5  of  Reference  [41],  Figure  3.15  shows  a  comparison  of  Lyttle’s  distribution  with 
that  computed  by  STABL.  A  difference  of  about  5  K  is  evident  along  the  frustum  of 
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Figure  3.15:  Surface  temperature  distribution  for  Stetson’s  blunt  cone  with  rn  =  3.81 
mm  as  computed  by  STABL  and  by  Lyttle  et  al.  [41]. 


the  cone.  The  reason  for  this  difference  is  currently  unknown,  but  it  could  be  related 
to  the  differences  in  thermal  conductivity  shown  in  Figure  3.5. 

A  solution  was  obtained  which  used  the  wall  temperature  distribution  computed 
by  Lyttle  as  a  fixed  boundary  condition,  rather  than  the  standard  adiabatic  wall 
condition.  The  adiabatic  temperature  distribution  computed  by  STABL  was  used 
forward  of  s  =  0.008  m  because  accurate  results  could  not  be  obtained  from  the  plot 
in  Reference  [41].  Figure  3.16  shows  the  computed  amplification  rates  using  each 
mean  flow  solution.  No  significant  differences  are  present  in  the  local  amplification. 
Figure  3.17  shows  the  N  factors  computed  with  each  wall  temperature  distribution. 
The  use  of  the  Lyttle  wall  temperature  distribution  results  in  a  slightly  lower  N 
factor  than  with  the  STABL  calculated  distribution,  which  is  opposite  the  expected 
trend  of  a  greater  second-mode  N  factor  for  a  colder  wall.  The  difference  appears  to 
be  in  the  critical  location,  as  the  local  slopes  of  the  two  curves  are  the  same. 
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Figure  3.16:  Amplification  rate  for  Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  = 
3.81  mm  computed  using  mean  flows  with  the  surface  temperature  distributions 
calculated  by  Lyttle  and  STABL. 


3.2.6  STABL  User  Options 

STABL  has  a  number  of  options  that  a  user  can  adjust  to  give  more  accurate 
results  or  to  match  the  methods  of  another  computation.  Many  combinations  of  var¬ 
ious  parameters  were  tested,  and  in  all  cases  the  effects  on  the  stability  results  were 
found  to  be  negligible.  Figure  3.18  and  Figure  3.19  show  the  effect  on  the  amplifica¬ 
tion  rate  and  N  factor  of  changing  two  of  these  options.  For  the  “without  chemistry” 
line,  chemical  reactions  and  vibrational  energy  modeling  were  turned  off.  This  is  not 
expected  to  have  a  significant  effect  on  the  results  for  wind  tunnel  conditions.  For 
the  “reduced  dissipation”  line,  the  numerical  dissipation  parameters  e *  and  e3  in 
DPLR  were  halved.  In  all  cases,  the  curves  overlay  each  other  completely,  showing 
the  independence  of  the  solution  with  respect  to  these  factors.  The  computational 
time  was  not  significantly  affected  by  either  of  these  changes. 
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Figure  3.17:  N  factors  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm  computed  us¬ 
ing  mean  flows  with  the  surface  temperature  distributions  calculated  by  Lyttle  and 
STABL. 


Figure  3.18:  Amplification  rate  curves  for  Stetson’s  blunt  cone  at  s/rn  =  175  with 
rn  =  3.81  mm  computed  using  various  user-specified  options  in  STABL. 
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Figure  3.19:  N  factor  curves  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm  computed 
using  various  user-specified  options  in  STABL. 


3.2.7  Freestream  Thermal  Non-Equilibrium 

All  of  the  other  computations  carried  out  in  this  chapter  assume  that  the  freestream 
flow  is  in  thermal  equilibrium.  This  is  consistent  with  the  other  computations  ref¬ 
erenced,  which  were  all  performed  with  codes  that  do  not  consider  any  effect  of  the 
vibrational  temperature. 

Roy  et  al.  [56]  showed  that  the  expansion  rate  in  the  Sandia  National  Laboratories 
hypersonic  wind  tunnel  with  the  Mach  8  nozzle  is  much  larger  than  the  thermal 
relaxation  rate,  causing  the  vibrational  temperature  of  the  freestream  flow  in  the 
test  section  to  be  frozen  near  the  stagnation  temperature.  It  seems  likely  that  a 
similar  situation  could  exist  in  AEDC’s  Tunnel  B,  which  has  a  similar  Mach  number 
and  slightly  higher  freestream  stagnation  temperature  than  the  633  K  found  in  the 
Sandia  tunnel.  Bertolotti  [57]  showed  that  vibrational  energy  relaxation  can  have  a 
significant  effect  on  stability.  He  states  that  since  the  vibrational  energy  does  not 
exceed  10%  of  the  total  internal  energy  of  air  until  the  flow  stagnation  temperature 
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Figure  3.20:  N  factors  calculated  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm  assum¬ 
ing  both  equilibrium  and  frozen  thermal  states. 


surpasses  800  K,  the  degree  of  thermal  nonequilibrium  does  not  significantly  impact 
stability  below  this  temperature.  Bertolotti  does  not  present  any  computations  for 
flows  with  stagnation  temperatures  below  1000  K. 

Since  stability  results  have  been  shown  to  be  very  sensitive  to  small  changes  in 
the  mean  flow,  it  seems  likely  that  there  could  be  some  effect  even  at  a  stagnation 
temperature  of  750  K.  The  degree  of  thermal  nonequilibrium  in  Tunnel  B  was  not 
calculated,  but  to  gauge  the  maximum  potential  effects  on  stability,  computations 
were  performed  with  the  freestream  vibrational  temperature  set  to  the  stagnation 
temperature.  The  true  vibrational  temperature  should  be  somewhere  between  the 
stagnation  temperature  and  the  freestream  static  temperature,  so  these  stability 
results  and  those  for  full  equilibrium  should  bracket  the  results  for  the  actual  condi¬ 
tions. 

Figure  3.20  shows  a  comparison  of  the  N  factors  obtained  when  thermal  equilib¬ 
rium  is  assumed  and  when  the  vibrational  temperature  is  frozen  at  the  stagnation 
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temperature.  Differences  are  evident  in  both  the  critical  location  for  a  given  fre¬ 
quency  and  the  growth  rate,  or  slope  of  the  curve.  When  the  thermal  state  is  frozen, 
N  =  5.5  is  first  reached  by  a  120  kHz  mode  (not  shown  on  figure)  at  s  =  0.959, 
which  is  4%  earlier  than  the  location  of  s  =  1.00  obtained  when  thermal  equilibrium 
is  assumed.  When  the  amplification  rate  at  s/rn  =  175  is  plotted  for  both  frozen  and 
equilibrium  flow,  no  difference  is  seen.  This  occurs  despite  the  difference  in  slope 
evident  in  Figure  3.20  at  both  150  and  170  kHz  at  that  location,  s  =  0.67  m.  This 
suggests  that  the  difference  in  slope  might  be  caused  by  changes  in  the  disturbance 
shape  factor,  which  are  included  in  a  PSE  analysis  but  not  in  a  local  LST  calculation. 

The  effect  on  the  amplification  rates  at  s/rn  =  175  is  many  orders  of  magnitude 
below  the  effect  on  the  N  factors.  The  reason  for  this  is  not  currently  known. 

3.2.8  Normalization 

All  of  the  references  cited  present  their  amplification  rates  in  non-dimensional 
form,  but  only  Lyttle  et  al.  [41]  state  the  numerical  values  used  in  the  normalization. 
The  others  all  use  the  length  scale  given  by  Equation  3.12. 


By  default,  PSE-Chem  will  output  both  dimensional  amplification  rates  and  those 
normalized  by  this  length  scale,  but  the  value  of  the  length  scale  depends  on  the  vis¬ 
cosity  model  employed.  Since  there  is  no  standard  viscosity  model  for  the  freestream 
temperature,  all  of  the  figures  in  this  thesis  are  presented  in  dimensional  quantities 
to  avoid  confusion.  The  data  digitized  from  the  references  was  subsequently  con¬ 
verted  to  dimensional  form  using  the  values  provided  for  the  Lyttle  et  al.  data  and 
Equation  3.12  for  all  of  the  others. 
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Figure  3.21:  The  effect  of  using  different  C  and  Fortran  compilers  on  complex  wave 
number  for  Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  =  3.81  mm. 


3.2.9  Compilers 

The  C  and  Fortran  compilers  used  to  compile  the  STABL  source  code  were  at 
one  time  found  to  make  a  difference  in  the  results  obtained.  The  documentation 
states  that  STABL  is  compatible  with  either  the  GNU  compilers  gcc  and  g77  or  the 
Portland  Group  compilers  pgcc  and  pgf  90  [26].  The  GNU  compilers  were  employed 
for  the  first  eight  months  the  code  was  used  because  of  their  ready  availability  on  a 
Linux  machine.  Differences  in  the  output  obtained  by  the  author  and  Heath  Johnson 
using  the  same  input  files  and  source  code  led  to  a  trial  of  the  Portland  Group 
compilers.  This  was  found  to  have  a  significant  effect,  as  illustrated  in  Figure  3.21, 
which  shows  the  different  results  obtained  with  each  set  of  compilers.  In  other  cases, 
using  the  GNU  compilers  caused  PSE-Chem  to  crash  without  producing  any  results 
at  all.  These  problems  were  not  encountered  with  the  Portland  Group  compilers. 

The  discrepancy  between  the  compilers  appears  to  have  been  caused  by  a  bug  in 
the  PSE-Chem  source  code  that  was  recognized  and  avoided  by  the  Portland  Group 


48 


Figure  3.22:  Amplification  curve  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm  at 
s  =  0.3  m  showing  the  existence  of  multiple  families  of  solutions  with  different  phase 
speeds. 


compilers  but  that  caused  failure  in  the  less  sophisticated  GNU  compilers.  That  bug 
has  since  been  fixed,  and  others  have  subsequently  used  the  GNU  compilers  without 
incident.  However,  this  is  an  example  of  the  importance  of  code  verification.  If 
the  same  conditions  had  not  been  run  by  multiple  users  on  multiple  machines  using- 
different  configurations,  the  bug  may  have  gone  unnoticed  indefinitely. 

3.2.10  Numerical  Behavior 

As  described  in  Section  2,  the  global  solution  procedure  produces  a  spectrum 
of  approximate  eigenvalues,  and  the  local  procedure  iterates  those  guesses  to  the 
converged  solution.  For  the  Stetson  et  al.  conditions,  and  presumably  others,  there 
exist  many  solutions  for  each  frequency,  and  different  guesses  will  converge  to  dif¬ 
ferent  solutions.  To  determine  the  complex  wavenumber  of  the  second  mode,  the 
local  solver  must  analyze  a  guess  that  converges  to  that  most  unstable  solution.  The 
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Figure  3.23:  Amplification  curves  for  Stetson’s  blunt  cone  at  s/rn  =  175  with  rn  = 
3.81  mm  obtained  with  the  original  extrapolation  routine  and  with  the  modified 
extrapolation  routine. 


solutions  can  be  differentiated  by  their  phase  speeds.  The  second  mode  typically  had 
phase  speeds  in  the  range  of  0.90-0.95,  whereas  the  other  solutions  would  have  phase 
speeds  greater  than  1.0.  Mack  shows  schematic  diagrams  of  the  different  families  of 
solutions  in  Figures  1  and  2  of  Reference  [58].  An  example  of  the  different  families 
of  results  for  the  Stetson  case  at  s/rn  =  175  is  shown  in  Figure  3.22. 

Mack  observes  that  the  radius  of  convergence  for  the  most  unstable  solution  can 
be  very  small  at  hypersonic  speeds,  and  that  was  observed  in  the  present  work.  As 
an  example,  PSE-Chem  gives  the  user  the  option  to  specify  a  wavenumber  guess  to 
the  local  solver.  Even  when  the  converged  results  of  a  previous  local  solution  at  the 
same  frequency  are  fed  back  to  the  solver  with  five  significant  figures,  often  the  new 
solution  would  converge  to  one  of  the  damped  families  of  solutions.  PSE-Chem  has 
options  to  extrapolate  guesses  for  new  frequencies  based  on  converged  solutions  at 
neighboring  frequencies.  Based  on  results  obtained  by  the  author,  the  extrapolation 
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Figure  3.24:  Spatial  growth  rate  and  N  factor  for  Stetson’s  blunt  cone  with  rn  =  3.81 
mm  computed  using  different  step  sizes  for  the  PSE  marching  procedure. 


algorithm  was  changed  to  enable  sweeps  of  both  increasing  and  decreasing  frequency, 
which  helped  to  ensure  that  once  the  second  mode  family  was  found  for  one  frequency 
it  would  be  found  for  the  rest  of  the  range.  Figure  3.23  shows  the  results  obtained 
when  the  original  extrapolation  routine  is  used  and  when  the  new  extrapolation 
routine  is  used. 

Another  numerical  difficulty  encountered  during  the  solution  process  was  an  in¬ 
stability  of  the  PSE  marching  procedure.  Although  the  pressure  gradient  is  modified 
in  the  PSE  marching  to  suppress  upstream  disturbances  [23],  the  solution  can  be¬ 
come  unstable  for  small  step  sizes.  In  PSE-Chem,  the  problem  manifested  itself  in 
oscillations  in  the  spatial  growth  rate,  as  shown  in  Figure  3.24.  PSE-Chem  provides 
an  option  to  the  user  to  skip  an  integer  number  of  axial  stations  during  each  march¬ 
ing  step.  This  allows  the  user  to  specify  a  large  number  of  axial  points  to  obtain 
a  more  accurate  mean  flow  without  sacrificing  the  numerical  stability  of  the  PSE 
procedure. 


51 


Figure  3.25:  Mean  velocity  profile  for  conditions  of  Stetson  et  al.  [33]  at  s/rn  =  175 
with  rn  =  3.81  mm  compared  with  the  results  of  Kufner  et  al.  [36]. 


3.3  Computational  Comparison 

The  mesh  used  for  the  final  computation  contained  440  axial  points  and  350  body- 
normal  points,  and  it  corresponds  to  the  440x350  (b)  line  in  Figure  3.14.  Chemical 
and  thermal  nonequilibrium  effects  were  modeled,  and  a  freestream  air  mixture  with 
standard  mass  fractions  of  76.7%  N2  and  23.3%  02  was  assumed.  The  freestream 
density  was  set  to  0.02464^-  to  match  the  conditions  thought  to  be  used  by  the  other 
researchers.  Thermal  equilibrium  was  assumed  in  the  freestream.  STABL  version 
1.26  was  used  for  both  the  mean  how  and  the  stability  calculations. 

Figure  3.25  shows  the  mean  tangential  how  velocity  at  s/rn  =  175.  Data  in  Figure 
la  of  Kufner  et  al.  [36]  was  digitized  using  the  method  described  in  Section  3.1  and 
is  plotted  on  Figure  3.25  for  comparison.  Agreement  is  generally  good,  although 
there  is  a  small  difference  near  the  edge  of  the  boundary-layer.  Figure  3.26  shows 
the  second  derivative  of  the  velocity  profile  with  respect  to  the  body-normal  distance 
at  the  same  station.  Data  from  Figure  lb  in  Reference  [36]  is  also  plotted.  Again, 
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Figure  3.26:  Second  normal  derivative  of  velocity  profile  for  conditions  of  Stetson  et 
al.  [33]  at  s/rn  =  175  with  rn  =  3.81  mm  compared  with  the  results  of  Kufner  et 
al.  [36]  and  Esfahanian  [35]. 


Figure  3.27:  Amplification  rate  for  conditions  of  Stetson  et  al.  [33]  at  s/rn  =  175 
with  rn  =  3.81  mm. 
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Table  3.4:  Location  of  the  peak  amplification  as  computed  by  STABL  and  several 
researchers. 


Computation 

co  (kHz) 

-cti  (1/m) 

STABL 

127 

16.5 

Malik  [34] 

132 

13.8 

Esfahanian  [35] 

129 

15.1 

Kufner  [36] 

126 

13.5 

Stilla  [36] 

122 

14.3 

Lyttle  [41] 

133 

14.7 

Zhong  [42] 

134 

15.2 

agreement  is  good  except  near  the  boundary-layer  edge.  Because  the  instabilities 
arise  near  the  boundary-layer  edge,  this  difference  could  be  significant. 

Figure  3.27  shows  the  amplification  curve  at  s/rn  =  175  compared  with  that  com¬ 
puted  by  other  researchers.  The  peak  frequency  is  slightly  lower  than  the  mean  peak 
frequency,  but  it  is  within  the  scatter  of  the  other  results.  The  peak  amplification 
rate  is  higher  than  any  of  the  other  researchers.  The  magnitude  of  the  difference  in 
the  peak  amplification  rate  is  comparable  to  the  scatter  in  the  other  results  shown. 
Table  3.4  summarizes  the  location  of  the  peaks  in  Figure  3.27. 

As  mentioned  previously,  most  of  the  literature  pertaining  to  the  Stetson  et  al. 
experiment  focuses  on  the  s/rn  =  175  axial  station.  Malik  et  al.  present  amplification 
rates  for  s/rn  =  215  in  Figure  15b  of  Reference  [34],  Figure  3.28  shows  a  comparison 
of  the  amplification  curve  computed  by  STABL  at  that  station  with  Malik’s  data. 
The  level  of  agreement  is  consistent  with  the  results  presented  in  Figure  3.27,  in 
that  the  peak  computed  by  PSE-Chem  is  shifted  up  and  to  the  left  with  respect  to 
Malik’s  peak. 
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Figure  3.28:  Amplification  rate  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm  at 
s/rn  =  215  compared  with  the  data  of  Malik  et  al.  [34], 


Figure  3.29:  Amplification  rate  curves  at  axial  locations  ranging  from  s  =  0.3  m  to 
s  =  1.0  m  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm. 
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Figure  3.30:  Maximum  N  factors  obtained  for  various  frequencies  for  Stetson’s  blunt 
cone  with  rn  =  3.81  mm. 


Figure  3.29  presents  the  local  amplification  rate  computed  at  axial  locations 
ranging  from  s  =  0.3  m  to  s  =  1.0  m.  No  amplification  is  found  forward  of  s  =  0.3 
m.  The  peak  amplification  grows  from  s  =  0.3  m  to  s  =  0.5  m  and  is  nearly  constant 
from  s  =  0.5  m  to  the  aft  end  of  the  cone.  The  most  amplified  frequency  decreases 
as  the  axial  distance  increases.  This  is  expected,  given  that  the  frequency  is  highly 
tuned  to  the  boundary-layer  thickness  [59].  One  interesting  phenomenon  is  that 
for  all  locations  from  s  =  0.4  m  aft,  amplification  occurs  for  frequencies  below  the 
primary  second  mode  frequency.  That  amplification  represents  a  2-D  first  mode 
instability,  and  amplification  in  that  frequency  range  is  greater  for  oblique  wave 
angles.  The  solution  from  70-150  kHz  at  s  =  0.3  m  represents  one  of  the  damped 
families  of  solutions  described  in  Section  3.2.10.  In  this  case,  that  solution  is  less 
stable  than  the  second  mode  solution  at  those  frequencies. 

Figure  3.30  displays  the  maximum  N  factors  obtained  for  various  frequencies 
for  Stetson’s  rn  =  3.81  mm  case.  To  obtain  these  N  factors,  a  test  matrix  was  con- 
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Table  3.5:  Transition  location  predicted  by  STABL  for  Stetson  et  al.  blunt  cone 
case  compared  with  that  predicted  by  other  researchers.  Transition  was  assumed  at 
N=5.5  for  all  computations. 


Computation 

Str  (m) 

Re 

J  toOO ,str 

STABL 

1.00 

8.192  x  106 

Esfahanian  and  Hejranfar  PNS  [40] 

1.00 

8.21  x  106 

Esfahanian  and  Hejranfar  IPNS  [40] 

0.973 

7.98  x  106 

Stilla  [37] 

0.956 

7.84  x  106 

Malik  et  al.  [34] 

0.956 

7.84  x  106 

Rosenboom  et  al.  [38] 

>  1.0 

>  8.2  x  106 

structed  consisting  of  all  combinations  of  the  starting  locations  s  =  0.1,  0.15,  ...,1.0m 
and  frequencies  uj  =  50,  60, ... ,  200  kHz.  PSE  marching  was  conducted  at  each  point 
in  the  test  matrix,  and  the  maximum  N  factors  for  each  frequency  are  plotted  in  Fig¬ 
ure  3.30.  If,  following  the  work  of  Malik  et  al.  [34],  Stilla  [37],  and  Esfahanian  and 
Hejranfar  [40],  transition  is  assumed  to  occur  at  N  =  5.5  in  this  noisy  environment, 
STABL  would  predict  transition  at  s  =  1.00  m  with  130  kHz  the  most  amplified  fre¬ 
quency.  Table  3.5  shows  a  comparison  of  the  transition  location  predicted  by  STABL 
with  that  of  the  other  researchers.  The  unit  Reynolds  number  specified  by  Stetson 
et  al.  of  2.50  x  106/ft  was  used  to  convert  between  the  Re and  str.  Agreement 
with  the  other  researchers  is  better  than  5%. 

Figure  3.31  shows  a  comparison  of  the  N  factors  computed  by  STABL  with  the 
local  N  factors  computed  by  Rosenboom  et  al.  [38].  Agreement  in  the  location  of 
amplification  of  the  various  frequencies  is  fair,  with  no  clear  trends  evident.  The 
STABL  140  kHz  and  160  kHz  lines  show  amplification  throughout  the  full  surface 
length,  in  contrast  to  the  Rosenboom  results,  which  begin  to  decay  at  s  =  0.85  and 
s  =  0.55  m,  respectively.  The  STABL  results  are  consistent  with  Figure  3.29,  which 
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Figure  3.31:  Comparison  of  N  factors  calculated  by  STABL  with  local  N  factors  of 
Rosenboom  et  al.  [38]  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm. 


shows  all  frequencies  up  to  180  kHz  amplifying  on  the  aft  portion  of  the  cone.  The 
reason  for  the  differences  is  not  clear.  Exact  results  for  Rosenboom  et  al.  are  not 
given  in  Table  3.5  because  the  data  do  not  extend  beyond  s  =  1.00  m  or  N  =  5.0, 
so  only  lower  bounds  for  the  predicted  transition  location  are  given. 

Figures  3.32  and  3.33  show  the  temperature  and  pressure  eigenfunctions  at  the 
aft  end  of  the  cone.  These  are  dimensional  quantities  obtained  for  an  initial  distur¬ 
bance  amplitude  of  0.001  K  or  Pa.  The  shape  of  the  curves  is  indicative  of  second 
mode  disturbances,  as  shown  by  Johnson  [23].  The  pressure  amplitude  shows  two 
peaks,  which  matches  the  results  of  Mack  [12]  for  second  mode  waves.  The  pressure 
amplitude  in  Figure  3.32  does  not  cross  zero  because  the  amplitude  is  plotted,  in 
contrast  to  Mack,  who  plotted  the  real  component  of  the  eigenfunction.  The  bump 
in  the  pressure  amplitude  and  phase  at  y/5  =  0.85  is  thought  to  be  a  numerical 
artifact,  but  it  was  not  investigated. 
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Figure  3.32:  Amplitude  and  phase  of  the  temperature  eigenfunction  at  the  aft  end 
of  the  cone  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm. 


Figure  3.33:  Amplitude  and  phase  of  the  pressure  eigenfunction  at  the  aft  end  of  the 
cone  for  Stetson’s  blunt  cone  with  rn  =  3.81  mm. 
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4.  First-Mode  Verification 

Two  cases  were  analyzed  to  verify  STABL’s  ability  to  correctly  calculate  first-mode 
disturbances.  Since  real  hypersonic  vehicles  have  some  degree  of  bluntness  to  control 
stagnation-point  heating,  a  blunt  cone  case  was  desired.  However,  no  previously 
calculated  blunt  cone  cases  for  which  the  first  mode  was  dominant  could  be  found. 
A  low  Mach  number  sharp  cone  was  found  that  had  been  shown  to  be  first  mode 
dominant.  In  addition,  a  large-bluntness  cone  was  found  for  which  the  first  mode 
had  a  degree  of  amplification  comparable  to  the  second  mode.  Both  cases  will  be 
discussed  in  this  chapter. 

4.1  Sharp  Cone  at  Mach  3.5 

4.1.1  Experimental  Conditions 

A  sharp  cone  was  analyzed  using  the  conditions  shown  in  Table  4.1  to  verify 
STABL’s  ability  to  calculate  first  mode  disturbances.  The  conditions  correspond 
to  experimental  data  obtained  by  Beckwith  et  al.  [60]  in  the  NASA  Langley  Mach 
3.5  pilot  quiet  tunnel.  This  tunnel  is  designed  to  give  freestream  disturbance  levels 
comparable  to  flight.  The  conditions  chosen  correspond  to  the  fourth  unit  Reynolds 
number  data  set  in  Run  5  of  Reference  [60].  The  first  mode  is  expected  to  be 
dominant  at  this  low  freestream  Mach  number  with  an  adiabatic  wall. 

This  case  was  previously  analyzed  by  Malik  as  case  QT1  in  Reference  [61]  using 
the  COSAL  stability  code.  This  older  LST-based  code  calculates  temporal  stability 
and  uses  Gaster’s  group  velocity  transformation  to  obtain  spatial  stability  results. 
The  mean  how  was  obtained  from  a  boundary-layer  code  and  provided  as  an  input 
to  COSAL. 
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Table  4.1:  Test  conditions  for  sharp  cone  at  Mach  3.5 


Condition 

Value 

Cone  half  angle 

5° 

Cone  length  (m) 

0.381 

Wall  temperature 

Adiabatic 

Fluid 

Air 

Moo 

3.5 

P O  (kPa) 

525 

T0  (K) 

319.0 

Too  (K) 

92.53 

Poo  (kg/m3) 

0.02592 

Reoo/m 

2.74xl07 

Sir  (m) 

0.278 

Roger  Kimmel  of  the  Air  Force  Research  Laboratory  (personal  communication, 
March-April  2005)  performed  additional  computations  for  this  case  using  the  eMahk 
stability  solver.  This  more  recent  code  is  also  widely  used  for  hypersonic  stability 
analysis.  Spatial  stability  results  are  calculated  directly,  and  a  perfect  gas  model  is 
assumed.  For  the  cases  shown,  the  mean  flow  was  provided  by  a  similarity  solver 
built  into  the  eMahk  code.  These  computations  were  performed  at  the  request  of  the 
author,  and  a  good  deal  of  collaboration  and  discussion  with  Dr.  Kimmel  occurred 
to  troubleshoot  differences  between  the  results. 

4.1.2  Computational  Comparison 


A  grid  with  450  axial  points  and  350  normal  points  was  used  for  the  DPLR2D 
mean  flow  solution.  Exponential  stretching  was  used  in  both  directions  to  cluster  grid 
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Figure  4.1:  Lines  of  constant  pressure  at  the  nosetip  for  a  sharp  cone  at  Mach  3.5. 


points  near  the  cone  tip.  This  clustering  increased  the  number  of  points  within  the 
boundary-layer  and  reduced  the  degree  of  numerical  error  caused  by  the  singularity 
at  the  nose.  Figure  4.1  shows  isobars  for  the  first  0.5%  of  the  cone  length.  The 
solution  is  valid  very  near  the  nose,  despite  the  singularity  caused  by  the  sharp  tip. 
The  number  of  grid  points  within  the  boundary-layer  varied  from  148  near  the  nose 
tip  to  109  at  the  base.  Chemical  and  thermal  nonequilibrium  effects  were  modeled 
with  thermal  equilibrium  assumed  in  the  freestream,  and  the  blended  viscosity  model 
was  used  for  all  calculations. 

Figure  4.2  shows  wall-tangential  mean  flow  velocity  profiles  computed  by  STABL 
and  eMahk  at  three  axial  locations.  The  velocities  computed  by  STABL  are  on  the 
order  of  10-20  m/s  larger  than  those  computated  by  eMahk  throughout  the  profile 
for  all  distances.  Similar  trends  are  seen  in  the  temperature  profiles,  shown  in  Fig¬ 
ure  4.3,  with  the  differences  being  on  the  order  of  5  K.  These  visible  differences  in  the 
temperature  and  velocity  can  be  expected  to  have  a  significant  effect  on  the  stability 
results.  Meaningful  second-derivative  comparisons  were  not  available  due  to  the  lack 
of  precision  in  the  output  of  the  eMahk  similarity  solution. 

Linear  stability  calculations  were  performed  at  several  axial  locations.  Figure  4.4 
shows  the  amplification  rate  plotted  as  a  function  of  both  frequency  and  wave  angle 
for  six  different  axial  locations.  The  contour  plots  shown  were  created  by  interpolat- 
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Figure  4.2:  Mean  flow  velocity  profiles  for  the  sharp  cone  at  Mach  3.5  at  several 
axial  locations. 


Figure  4.3:  Mean  flow  temperature  profiles  for  the  sharp  cone  at  Mach  3.5  at  several 
axial  locations. 
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ing  between  the  results  of  twelve  solutions  at  each  location.  Each  of  the  individual 
solutions  had  a  varying  frequency  but  constant  /3.  The  constant  dimensional  span- 
wise  wavenumber  condition  obeys  the  irrotationality  condition  on  the  wavenumber 
vector  described  by  Mack  [12,58].  The  wave  angle  was  computed  for  every  point 
using  Equation  4.1. 

ip  =  arctan  ^  (4-1) 

Several  qualitative  conclusions  can  be  drawn  from  Figure  4.4.  The  most  amplified 
frequency  decreases  as  R  increases.  The  maximum  amplification  rate  also  decreases 
as  R  increases,  although  that  effect  is  not  seen  very  near  the  nose.  The  wave  angle 
of  the  most  amplified  mode  increases  from  approximately  60°  very  near  the  nose  to 
approximately  70°  on  the  aft  portion  of  the  cone. 

Several  mean  flow  solutions  using  different  numbers  of  grid  points  within  the 
boundary-layer  were  analyzed  to  ensure  grid-independence  of  the  solution.  Figure  4.5 
shows  the  effect  of  the  different  grids  on  the  N  factor  calculations.  These  calculations 
are  for  a  mode  at  w  =  104  kHz  and  (3  =  2400/m  with  marching  beginning  at 
s  =  0.0339  m.  Four  of  the  solutions  are  very  similiar;  the  only  difference  is  seen 
for  the  grid  with  72  points  in  the  boundary-layer.  The  reason  for  this  repeatable 
difference  is  not  clear.  The  grid  with  109  points  in  the  boundary-layer  was  used  for 
all  other  calculations  in  this  section. 


To  obtain  the  maximum  N  factors  for  this  case,  a  test  matrix  was  constructed 
consisting  of  u  =  25,  50, ... ,  250  kHz,  [3  =  1000,  2000, . . . ,  8000  1/m,  and  starting 
location  s  =  0.02,  0.04, . . . ,  0.12  m.  The  modes  producing  the  maximum  N  factors 
at  each  location  are  shown  in  Figure  4.6.  Amplification  began  at  s  =  0.02  m  for  the 
higher  frequencies  and  shortly  thereafter  for  the  75  kHz  mode.  At  the  experimental 
transition  location  of  s  =  0.278  m,  the  most  amplified  mode  has  a  frequency  of  75 
kHz,  a  (3  of  2000  1/m,  a  wave  angle  of  64°,  and  a  N  factor  of  12.2.  This  N  factor 
is  larger  than  the  range  of  9-11  commonly  calculated  at  experimentally  measured 
transition  locations  in  low  disturbance  environments.  For  comparison,  the  value  of 
Reg/Me  is  154  at  the  transition  location,  as  shown  in  Figure  4.7. 


(e)  R=2000  (s  =  0.136  m) 


Figure  4.4:  Amplification  rate  for  the 
frequency  and  wave  angle  for  various  ax 
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Figure  4.5:  N  factor  calculations  for  the  sharp  cone  at  Mach  3.5  with  five  separate 
mean  flow  grids.  Calculations  are  for  a  mode  with  co  =  104  kHz  and  /3  =  2400/m, 
beginning  at  s  =  0.0339  m. 


Figure  4.6:  Maximum  N  factors  obtained  from  many  combinations  of  frequency, 
spanwise  wavenumber,  and  starting  location  for  the  sharp  cone  at  Mach  3.5. 


Figure  4.7:  Ratio  of  edge  momentum  thickness  Reynolds  number  to  edge  Mach 
number  for  the  sharp  cone  at  Mach  3.5. 


Figure  4.8:  The  phase  speed  for  the  five  most  amplified  modes  on  the  sharp  cone  at 
Mach  3.5. 
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Figure  4.9:  The  wave  angle  for  the  five  most  amplified  modes  on  the  sharp  cone  at 
Mach  3.5. 


Figure  4.8  shows  the  phase  speed  for  each  of  the  five  modes  shown  in  Figure  4.6. 
The  phase  speed  is  always  within  the  range  0.55-0.85,  which  is  lower  than  the  0.90- 
1.00  range  found  for  second  mode  waves.  The  phase  speed  increases  for  each  mode 
as  the  surface  distance  increases,  as  is  expected  based  on  Equation  3.11  since  ar 
decreases  as  the  boundary-layer  thickness  increases.  The  phase  speed  increases  at 
all  distances  as  the  frequency  increases,  as  is  also  expected  based  on  Equation  3.11. 
The  phase  speed  of  the  most  amplified  mode  at  the  transition  location  is  0.73. 

Figure  4.9  shows  the  wave  angles  for  each  of  the  five  modes  shown  in  Figure  4.6. 
The  wave  angle  increases  as  the  surface  distance  increases,  which  is  expected  based 
on  Equation  4.1  since  ar  decreases  and  j3  is  held  constant.  No  clear  trend  of  wave 
angle  with  frequency  is  evident.  The  wave  angle  of  the  most  amplified  mode  at  the 
transition  location  is  64°. 

Figure  4.10  shows  the  amplitude  and  phase  of  the  temperature  eigenfunction  at 
the  aft  end  of  the  cone  for  a  104  kHz  disturbance  with  (3  =  2400  1/m.  Figure  4.11 
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Figure  4.10:  Amplitude  and  phase  of  the  temperature  eigenfunction  for  a  104  kHz 
disturbance  with  j3  =  2400  1/m  at  the  aft  end  of  the  sharp  cone  at  Mach  3.5. 


Figure  4.11:  Amplitude  and  phase  of  the  pressure  eigenfunction  for  a  104  kHz  dis¬ 
turbance  with  j3  =  2400  1/m  at  the  aft  end  of  the  sharp  cone  at  Mach  3.5. 
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shows  the  pressure  eigenfunction  for  the  same  disturbance.  The  one  peak  present  in 
both  the  temperature  and  pressure  amplitudes  is  indicative  of  a  first  mode  distur¬ 
bance.  The  shape  of  the  pressure  amplitude  curve  matches  well  with  that  shown  by 
Mack  [12]  for  first  mode  waves.  The  amplitudes  shown  are  for  an  initial  disturbance 
amplitude  of  0.001  at  s  =  0.0339  m. 

4.1.3  Verification  and  Validation  Issues 

Comparisons  were  made  between  the  present  results  and  those  obtained  by  Malik 
using  COSAL  [61]  and  Kimmel  using  eMahk  (personal  communication,  March-April 
2005).  Malik  calculated  N  factor  growth  rates  for  six  frequencies  ranging  from  22-188 
kHz.  His  most  amplified  mode  at  the  experimental  transition  location  has  a  frequency 
of  78  kHz,  a  wave  angle  of  approximately  60°,  and  a  N  factor  of  10.1.  Compared  to 
the  transition  N  factor  calculated  by  STABL  of  12.2,  this  is  a  significant  difference, 
although  the  frequency  and  wave  angle  of  the  most  amplified  disturbances  agree  to 
within  10%. 

Figure  4.12  shows  a  comparison  of  the  N  factors  calculated  by  STABL  and  those 
calculated  by  Kimmel  using  the  eMahk  code  for  a  mode  at  104  kHz  with  j3  =  2400 
1/m.  The  N  factors  calculated  by  STABL  are  approximately  30%  higher  than  those 
calculated  by  eMahk.  The  agreement  is  much  worse  than  that  seen  in  Figures  4.2 
and  4.3  for  the  mean-flow  profiles.  The  difference  in  N  factor  may  be  due  to  the 
subtle  differences  in  the  mean  flow,  or  it  may  be  due  to  differences  in  the  stability 
solvers.  An  attempt  was  made  to  analyze  the  eMahk  meanflow  using  PSE-Chem, 
but  publication-quality  results  were  not  obtained  due  to  problems  converting  the 
similarity  solution  to  a  general  2D  flowfield. 

4.2  Large-Bluntness  Cone  at  Mach  8 

All  hypersonic  vehicles  have  some  degree  of  nosetip  bluntness  to  control  stagna¬ 
tion  point  heating.  For  a  blunt  cone,  when  the  nose  radius  increases  for  a  fixed,  high 
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Figure  4.12:  N  factors  calculated  using  STABL  compared  with  results  of  Kimmel 
using  eMahk  for  a  mode  at  104  kHz  with  j3  =  2400  1/m  for  the  sharp  cone  at  Mach 
3.5. 


freestream  Mach  number,  the  first  mode  may  become  dominant.  For  this  reason, 
a  very  blunt  cone  was  analyzed  using  the  conditions  shown  in  Table  4.2  to  further 
verify  first  and  second-mode  computations.  The  flow  conditions  match  those  of  the 
Stetson  et  al.  [33]  blunt  cone  experiment,  but  the  bluntness  was  increased  while 
the  ratio  of  the  nose-tip  radius  to  the  body  length  was  kept  constant,  making  this 
case  impractical  for  wind  tunnel  experiments.  This  case  was  previously  analyzed  by 
Rosenboom  et  al.  [38,62],  This  case  was  chosen  for  first  mode  verification  because 
Rosenboom  et  al.  found  first  mode  N  factor  growth  at  a  level  comparable  to  second 
mode  growth.  Their  LST  based  calculations  showed  a  maximum  first  mode  N  factor 
at  s  =  11  m  of  approximately  6.75  for  a  2  kHz  disturbance.  For  comparison,  the 
maximum  second  mode  N  factor  Rosenboom  et  al.  calculated  was  8.6  at  34  kHz. 
Kufner  [63]  performed  calculations  for  this  case  without  finding  any  first  mode  in¬ 
stabilities.  However,  Rosenboom  points  out  that  his  second  mode  calculations  for 
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Table  4.2:  Test  conditions  for  cone  with  rn  =  42.67  mm  at  Mach  8 


Condition 

Value 

Nose  radius  (mm) 

42.67  mm 

Cone  half  angle 

7° 

Maximum  arc  length  (m) 

11.35 

Wall  temperature 

Adiabatic 

Fluid 

Air 

Moo 

7.99 

P O  (kPa) 

4000 

T0  (K) 

750 

Too  (K) 

54 

Poo  (kg/m3) 

2.7xl0-2 

Reoo/m 

2.5  xlO6 

this  case  suffered  from  numerical  difficulties  that  manifested  themselves  as  wiggles 
in  the  disturbance  growth  rates.  This  casts  doubt  on  the  first  mode  calculations  for 
the  present  case. 

An  accurate  mean  flow  was  computed  using  a  mesh  with  450  axial  points  and 
350  wall-normal  points.  Approximately  100  grid  points  were  contained  within  the 
boundary-layer  over  the  cone  frustum.  The  edge  Mach  number  was  4  near  the  cone 
shoulder  and  increased  to  6.3  far  downstream  of  the  nose. 

4.2.1  First  Mode  Stability 

Rosenboom  et  al.  [38]  calculated  the  first  mode  critical  point  to  be  0.15  m.  Figure 
23  of  Reference  [38]  shows  all  frequencies  from  1-9  kHz  becoming  unstable  around 
0.15  m,  followed  by  a  region  of  approximately  5  meters  where  they  again  become 
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Figure  4.13:  First  mode  N  factors  for  the  large  bluntness  cone  at  Mach  8. 


stable,  and  a  region  from  s  =  7  m  to  the  end  where  all  frequencies  from  1-9  kHz 
are  again  unstable.  Figure  A.  12  of  Reference  [62]  shows  the  wave  angles  for  the 
instabilities  Rosenboom  calculated.  For  frequencies  from  1-5  kHz,  the  wave  angle  is 
greater  than  80°  at  the  critical  point,  drops  to  around  60°  almost  immediately,  and 
eventually  levels  off  between  65°  and  70°.  For  frequencies  greater  than  5  kHz,  the 
wave  angle  is  less  than  50°  at  the  critical  point  but  rises  to  the  trend  of  the  lower 
frequencies  by  s  =  1  m.  Figure  28  of  Reference  [38]  shows  N  factors  for  this  case. 
Rosenboom  shows  amplification  from  0. 5-2.0  m,  a  short  period  of  damping,  and  then 
monotonic  growth  from  approximately  4  m  to  the  end  of  the  cone.  The  figure  was 
digitized  using  the  Un-Scan-It  software,  and  N  =  5.5  is  first  reached  at  s  =  9.7  m 
by  the  2  kHz  mode  at  a  wave  angle  of  70°. 

Computations  were  performed  with  STABL  to  compare  with  the  results  of  Rosen¬ 
boom  et  al.  Figure  4.13  shows  samples  of  the  first  mode  N  factors  calculated 
by  STABL.  A  test  matrix  consisting  of  combinations  of  s  =  0.5, 1.0, ... ,  10.5  and 
j3  =  20, 40, . . . ,  100  1/m  was  analyzed.  Combinations  involving  additional  values 


73 


Table  4.3:  Transition  location  predictions  for  the  large  bluntness  cone  at  Mach  8. 
Transition  is  assumed  at  N  =  5.5. 


Computation 

Sir  (m) 

io  (kHz) 

First  Mode 

STABL 

9.7 

3.33 

o 

O 

Rosenboom  et  al.  [38,62] 

9.7 

2 

o 

O 

Second  Mode 

STABL 

9.0 

34 

0° 

Rosenboom  et  al.  (local)  [38] 

9.6 

34 

0° 

Rosenboom  et  al.  (nonlocal)  [38] 

9.8 

34 

0° 

of  j3  <  20  1/m  were  also  analyzed  separately,  but  no  significant  amplification  was 
found.  For  all  cases  shown,  100  frequencies  within  the  range  1-10  kHz  were  analyzed 
at  the  starting  location  using  LST,  and  the  lowest  unstable  frequency  was  used  for 
PSE  marching.  For  these  particular  conditions,  that  lowest  unstable  frequency  was 
always  the  lowest  for  which  the  global  and  local  solvers  could  obtain  converged  solu¬ 
tions.  The  spanwise  wavenumber  (3  and  the  frequency  u  were  fixed  as  the  marching 
progressed  downstream.  Table  4.3  summarizes  the  comparison.  N  =  5.5  is  first 
reached  at  s  =  9.7  m  by  a  3.33  kHz  instability  at  70°.  Despite  differences  in  the  fre¬ 
quency,  the  transition  location  and  wave  angle  predicted  by  STABL  and  Rosenboom 
et  al.  are  the  same.  In  addition,  the  general  characteristics  of  the  N  factor  curves 
calculated  by  STABL  and  Rosenboom  et  al.  agree  well.  Both  show  an  initial  hump, 
followed  by  brief  damping  and  a  long  rise. 
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Figure  4.14:  Second  mode  amplification  curves  for  the  large-bluntness  cone  at 
Mach  8. 


4.2.2  Second  Mode  Stability 

A  comparison  was  also  made  with  the  second  mode  stability  results  of  Rosenboom 
et  al.  Rosenboom  shows  second  mode  amplification  beginning  at  s  =  5.5  m  at  40 
kHz.  N  =  5.5  is  reached  at  approximately  9.6  m  for  the  local  calculation  and 
approximately  9.8  m  for  the  nonlocal  calculation  at  a  frequency  of  34  kHz  for  both 
local  and  nonlocal. 

Figure  4.14  shows  the  second  mode  amplification  curves  calculated  by  STABL. 
For  the  computations  shown,  a  frequency  range  of  20-60  kHz  was  analyzed,  but  no 
eigenvalues  were  found  by  the  global  solver  at  frequencies  greater  than  35  kHz  for  any 
location.  These  calculations  found  a  critical  frequency  of  34  kHz  at  s  =  6.5  m.  The 
peak  amplification  increases  and  the  unstable  frequencies  decrease  as  s  increases.  The 
critical  location  is  slightly  aft  of  Rosenboom’s  critical  location  of  5.5  m  at  40  kHz. 
However,  since  the  lack  of  eigenvalues  above  35  kHz  in  PSE-Chem’s  calculations 
appears  to  have  been  caused  by  unknown  numerical  difficulties,  higher  frequency 
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Figure  4.15:  Second  mode  N  factors  for  the  large-bluntness  cone  at  Mach  8. 


modes  may  have  provided  better  agreeement  in  the  critical  location.  Rosenboom 
shows  the  34  kHz  mode  beginning  to  amplify  at  s  =  7  m,  which  is  comparable  to 
the  location  calculated  by  STABL  of  6.5  m. 

Figure  4.15  shows  the  second  mode  N  factor  curves  generated  by  STABL.  A 
frequency  range  of  20-60  kHz  was  analyzed  in  2  kHz  increments  for  the  global  solution 
and  0.2  kHz  increments  for  the  local  solution,  but  no  eigenvalues  were  found  by  the 
global  solver  at  frequencies  greater  than  35  kHz  for  any  location.  STABL  selected 
the  critical  frequency  at  each  starting  location  to  begin  PSE  marching.  Table  4.3 
summarizes  the  comparison.  Amplification  begins  at  s  =  6.5  m,  slightly  behind 
the  critical  location  of  s  =  5.5  m  calculated  by  Rosenboom  et  al.  An  N  factor  of 
5.5  is  first  reached  at  s  =  9.0  m  by  a  34  kHz  wave.  This  is  slightly  ahead  of  the 
location  calculated  by  Rosenboom  of  s  =  9.6  m  for  the  local  calculation  or  s  =  9.8 
for  the  nonlocal  calculation,  but  the  difference  is  on  the  order  of  the  starting  distance 
resolution  of  0.5  m.  Additional  computations  should  be  performed  to  reduce  the  gap 
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Figure  4.16:  A  comparison  of  the  spatial  growth  rates  of  the  largest  first  and  second¬ 
mode  N  factors  for  the  large-bluntness  cone  at  Mach  8. 


between  STABL  computations.  Overall  second-mode  agreement  between  STABL 
and  Rosenboom  is  good. 

An  outstanding  question  in  stability  theory  is  the  nature  of  the  interaction  be¬ 
tween  the  first  and  second  modes.  Figure  4.16  shows  the  amplification  rates  of  the 
most  amplified  first  and  second  mode  waves,  and  Figure  4.17  shows  the  N  factors  for 
those  waves.  Both  figures  show  the  edge  Mach  number  distribution.  The  apparent 
discontinuity  in  the  edge  Mach  number  at  s  =  3  m  is  caused  by  numerical  difficulties 
of  the  type  described  in  Section  3.2.3.  The  first  mode  amplifies  over  the  full  length 
of  the  body.  The  second  mode  begins  to  amplify  when  the  edge  Mach  number  is  5.7, 
and  it  becomes  dominant  when  the  edge  Mach  number  is  6.0. 

The  calculations  also  showed  another  unstable  mode  that  was  not  expected.  Fig¬ 
ure  4.18  shows  the  N  factors  calculated  by  two-dimensional  waves  at  1  kHz.  The 
different  curves  represent  different  starting  locations.  The  waves  amplify  very  slightly 
until  s  =  3  m,  at  which  point  they  amplify  very  rapidly  until  s  =  4  m.  They  then 


Figure  4.17:  A  comparison  of  the  largest  first  and  second-mode  N  factors  calculated 
by  STABL  for  the  large  bluntness  cone  at  Mach  8. 


s  (m) 

Figure  4.18:  N  factors  of  a  two-dimensional  instability  at  1  kHz  on  the  large-bluntness 
cone  at  Mach  8.  The  four  lines  represent  different  beginning  marching  locations. 
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decay  slightly  until  s  =  10  m,  at  which  point  they  again  amplify  rapidly  until  the 
end  of  the  body.  N  =  5.5  is  first  reached  at  s  =  10.6  m.  These  waves  were  found 
serendipitously  during  a  search  over  the  1-10  kHz  frequency  range  when  the  lowest 
amplified  frequency  was  chosen  for  marching.  Although  this  search  was  run  every 
0.5  m,  only  the  starting  locations  between  1.5  and  3.0  m  exhibited  this  behavior. 
The  reason  for  the  different  amplification  between  the  different  starting  locations  is 
unknown.  This  might  be  a  new  phenomenon,  but  further  investigation  is  necessary 
before  it  can  be  characterized. 
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5.  Comparisons  to  Purdue  Mach-6  Ludwieg  Tube 

Experiments 

As  discussed  in  Chapter  3,  the  Stetson  et  al.  [33]  hot-wire  measurements  are  of  lim¬ 
ited  utility  for  validation  purposes.  In  addition,  the  expense  of  conducting  research 
experiments  in  a  production  wind  tunnel  like  Tunnel  B  makes  it  unlikely  that  the  test 
can  be  replicated  in  that  tunnel.  For  this  reason,  Rufer  and  Schneider  [64]  have  been 
developing  equipment  and  procedures  to  conduct  hot-wire  measurements  of  instabil¬ 
ity  waves  in  the  Purdue  Mach-6  Ludwieg  Tube  (M6LT).  Although  this  wind  tunnel  is 
designed  for  quiet  flow,  it  is  currently  quiet  only  at  low  unit  Reynolds  numbers  [11], 
and  the  experiments  described  here  take  place  under  conventional-noise  conditions. 
The  measurement  of  natural  second-mode  instability  wave  growth  is  complicated 
by  the  six-to-ten-second  usable  run  time  of  the  M6LT,  which  makes  calibration  in 
the  tunnel  difficult  and  results  in  frequent  wire  breakage  from  the  many  start-up 
and  shut-down  cycles.  This  chapter  will  present  computations  that  are  designed  to 
match  the  conditions  of  several  of  Rufer’s  experiments. 

The  test  conditions  for  the  cases  discussed  here  are  summarized  in  Table  5.1. 
The  total  pressure  corresponds  to  the  pressure  in  the  tunnel  driver  tube,  and  it  was 
set  to  one  of  the  three  values  shown  in  the  table.  The  specific  value  will  be  noted 
for  each  case.  In  addition,  as  the  driver  tube  empties,  the  total  pressure  decreases 
throughout  the  run,  and  the  value  given  in  Table  5.1  is  the  nominal  total  pressure  at 
the  beginning  of  the  run.  As  the  hot-wire  measurements  were  taken  within  the  first 
second  of  operation,  the  actual  pressure  should  be  close  to  the  nominal  pressure.  The 
freestream  Mach  number  also  changes  as  the  total  pressure  changes  throughout  the 
run.  The  wall  temperature  listed  in  Table  5.1  is  an  estimate.  Although  the  driver 
tube  section  of  the  M6LT  is  heated  to  433  K,  the  temperature  of  the  test  section  is 
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Table  5.1:  Nominal  test  conditions  for  the  hot-wire  experiments  of  Rufer  [64] 


Condition 

Sharp  Cone 

Blunt  Cone 

Half  angle 

7° 

7° 

Length,  x  (m) 

0.5689 

0.5653 

Base  diameter  (m) 

0.1397 

0.1397 

Nose  radius  (mm) 

0 

0.508 

Wall  temperature  (K) 

300 

300 

Fluid 

Air 

Air 

Moo 

5.8 

5.8 

P O  (psia) 

90,  125 

45,  125 

T0  (K) 

433 

433 

Too  (K) 

56 

56 

CO 

to 

8 

0.03008,  0.04177 

0.01504,  0.04177 

believed  to  be  at  or  near  room  temperature.  The  short  run  time  is  assumed  to  be 
insufficient  to  cause  significant  heating  of  the  model.  The  sensitivity  of  the  stability 
to  small  deviations  from  these  assumptions  will  be  discussed. 

Distances  in  this  section  are  expressed  using  both  the  axial  length  x  and  the  arc 
length  s.  STABL  uses  the  arc  length  for  all  of  its  calculations,  but  the  experimentalist 
uses  x.  These  are  related  by  Equation  5.1,  where  the  cone  half  angle  4>  is  given  in 
degrees. 


sW  =  2*r„^  +  1  “  r”(1  “ 


for  all  x  >  rT, 


(5.1) 


360  cos  < 

For  a  sharp  cone,  rn  =  0  and  Equation  5.1  reduces  to  s(x)  =  xj  cos  <f). 

Mean  flows  were  computed  for  all  cases  using  grids  with  440  axial  points  and  350 
streamwise  points.  The  outer  grid  boundary  was  carefully  specified  to  closely  match 
the  angle  of  the  shock.  Stretching  was  used  to  cluster  the  points  near  the  nose  and 
the  body  surface.  Approximately  100  wall-normal  points  were  contained  within  the 
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Figure  5.1:  N  factors  calculated  for  Rufer’s  blunt  cone  at  125  psia  using  two  separate 
grids. 


boundary  layer.  A  second  mean  flow  was  calculated  using  a  440x300  mesh  for  one 
of  the  cases.  Although  the  outer  grid  boundary  was  the  same  for  each  grid,  different 
stretching  parameters  were  used.  The  second  grid  had  a  miminum  of  70  points  in 
the  boundary  layer.  Figure  5.1  shows  three  sets  of  N  factor  calculations  using  each 
of  the  grids.  The  results  are  almost  identical,  which  suggests  that  the  results  are 
grid  independent.  The  results  of  this  test  and  those  shown  in  previous  chapters  were 
deemed  sufficient  to  assume  grid  independence  for  the  other  cases  calculated  in  this 
chapter. 

5.1  Sharp  Cones 

5.1.1  Sharp  Cone:  90  psia 

Figure  5.2  shows  a  comparison  of  the  boundary-layer  profiles  calculated  by  STABL 
and  measured  by  Rufer  at  two  axial  locations.  Calibrated  hot-wire  measurements 
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Figure  5.2:  A  comparison  of  the  mass  flux  calculated  by  STABL  with  the  uncal¬ 
ibrated  hot  wire  measurements  of  Rufer  at  two  locations  on  the  sharp  cone  at  90 
psia. 


are  needed  for  a  direct  comparison  between  experiment  and  computation.  As  of 
this  writing,  only  uncalibrated  hot-wire  measurements  have  been  obtained.  As  the 
hot-wire  mean  voltage  is  proportional  to  the  mass  flux,  a  limited  comparison  can 
still  be  made  using  these  uncalibrated  measurements.  The  height  and  shape  of  the 
boundary-layer  profiles  can  be  compared  by  adjusting  the  scales  of  the  computa¬ 
tions  and  experiments  to  make  the  curves  fit  as  well  as  possible.  However,  a  full 
comparison  will  require  calibrated  measurements.  The  height  of  the  hot  wire  rela¬ 
tive  to  the  cone  wall  is  determined  using  a  telescope  and  custom-built  lens  system 
that  is  designed  to  correct  for  the  optical  distortion  caused  by  the  thick,  curved 
plexiglass  window.  The  bias  error  in  the  measured  heights  is  estimated  to  be  within 
0.15  mm.  The  distance  between  successive  data  points  is  controlled  by  the  traverse 
mechanism,  and  the  error  is  presumed  to  be  negligible.  Measurements  were  taken 
on  separate  runs  to  minimize  the  effect  of  the  pressure  drop  during  the  run.  For  this 
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Figure  5.3:  Amplification  rate  curves  for  Rufer’s  sharp  cone  at  90  psia.  The  spacing 
between  the  curves  is  two  inches. 


case,  the  boundary  layer  thickness  calculated  by  STABL  and  measured  by  Rufer  ap¬ 
pear  to  agree  to  within  the  uncertainty  in  the  measurements  for  both  axial  locations. 
Agreement  with  the  profile  shape  is  good,  particularly  for  the  x  =  11.75  inch  data. 
Separate  scales  are  used  for  the  voltage  for  each  axial  distance,  but  for  clarity  only 
the  scale  for  the  x  =  11.75  inch  measurements  is  shown. 

Figure  5.3  shows  amplification  rate  curves  calculated  by  STABL  for  several  axial 
locations.  As  expected,  the  most  unstable  frequency  decreases  as  the  boundary 
layer  thickens.  The  peak  amplification  rate  also  decreases  with  increasing  distance 
downstream.  This  trend  is  similar  to  that  seen  in  Figures  3.29  and  4.4. 

N  factor  curves  for  this  case  are  shown  in  Figure  5.4.  Modes  with  frequencies 
higher  than  300  kHz  are  not  shown,  as  they  were  not  necessary  to  define  the  N-factor 
envelope  in  the  region  near  transition.  Figure  14  in  Rufer  and  Schneider  [64]  shows 
suspected  second-mode  instabilities  at  x  =  11.75  inches  and  x  =  14.0  inches  but  not 
at  x  =  16.7  inches.  This  is  taken  as  an  indication  that  transition  occurred  somewhere 
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Figure  5.4:  N  factors  for  Rufer’s  sharp  cone  at  90  psia.  The  spacing  between  the 
curves  is  10  kHz. 


between  14.0  and  16.7  inches,  although  it  should  be  noted  that  this  may  not  coincide 
with  the  transition  location  that  would  be  found  using  a  different  detection  method. 
Using  the  envelope  of  the  curves  shown  in  Figure  5.4,  transition  occurred  between  N 
factors  of  4.3  and  4.7.  The  expected  transition  N  factor  will  be  lower  for  experiments 
in  the  M6LT  than  for  those  in  Tunnel  B  because  the  noise  level  is  higher  in  smaller 
tunnels  [10].  Thus,  transition  N  factors  less  than  5  are  not  surprising.  The  most 
amplified  frequencies  calculated  by  STABL  at  x  =  11.75  inches  and  x  =  14.0  inches 
are  220  kHz  and  200  kHz.  These  are  within  10%  of  the  peaks  on  the  amplification 
curve  shown  in  Figure  14  of  Rufer  and  Schneider  [64],  which  are  235  kHz  and  200 
kHz,  respectively.  Table  5.2  summarizes  the  comparison  with  Rufer’s  data  for  all  of 
the  cases  discussed  in  this  chapter. 
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Table  5.2:  Comparison  of  most  unstable  frequencies  calculated  by  STABL  and  mea¬ 
sured  by  Rufer  [64],  Transition  N  factors  are  also  shown  for  each  case. 


Case 

Location  (in) 

Rufer  uj  (kHz) 

STABL  u  (kHz) 

Ntr 

Sharp  90 

11.75 

235 

220 

4. 3-4. 7 

14.0 

200 

200 

Sharp  125 

10.3 

280 

220 

>  5.1 

12.3 

250 

280 

Blunt  125 

14.25 

240 

260 

>  5.4 

Blunt  45 

18.0 

150 

140 

3. 6-4.0 

20.0 

130 

130 

Figure  5.5:  Amplification  rate  curves  for  Rufer’s  sharp  cone  at  125  psia.  The  spacing 
between  the  curves  is  two  inches. 
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Figure  5.6:  N  factors  for  Rufer’s  sharp  cone  at  90  psia.  The  spacing  between  the 
curves  is  20  kHz. 


5.1.2  Sharp  Cone:  125  psia 

Figures  10  and  11  of  Reference  [64]  show  hot-wire  spectra  on  a  sharp  cone  with 
an  initial  driver-tube  pressure  of  125  psia.  Although  the  cone  used  in  this  experiment 
was  shorter  than  that  used  in  the  others,  the  same  computational  grid  was  used  as 
for  the  90  psia  sharp  cone.  Figure  5.5  shows  the  amplification  rate  curves  calculated 
by  STABL  for  this  case,  and  Figure  5.6  shows  the  N  factors.  The  computation  was 
carried  out  with  10  kHz  spacing,  but  only  20  kHz  spacing  is  shown  in  the  figure  for 
clarity.  Rufer’s  frequency  spectra  show  a  peak  amplitude  of  approximately  280  kHz 
at  x  =  10.3  inches  and  250  kHz  at  x  =  12.3  inches.  STABL  calculates  an  N  factor 
of  4.5  at  300  kHz  at  x  =  10.3  inches,  and  an  N  of  5.1  at  280  kHz  at  x  =  12.3  inches. 
Transition  was  not  observed  for  this  condition,  so  the  transition  N  factor  is  presumed 
to  be  greater  than  5.1. 
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Figure  5.7:  A  comparison  of  the  mass  flux  calculated  by  STABL  with  the  uncali¬ 
brated  hot  wire  measurements  of  Rufer  at  x  =  14.25  inches  on  the  blunt  cone  at  125 
psia. 


5.2  Blunt  Cones 

5.2.1  Blunt  Cone:  125  psia 

Figure  5.7  shows  a  comparison  of  the  boundary-layer  profiles  measured  by  Rufer 
and  computed  by  STABL  for  the  blunt  cone  with  a  125  psia  driver  tube  pressure. 
These  measurements  were  conducted  before  the  telescope  system  was  in  place  to 
increase  the  accuracy  of  the  hot-wire  positioning  relative  to  the  cone  wall.  Thus, 
although  the  distance  between  Rufer’s  points  is  accurate,  the  entire  set  may  have  a 
constant  error  in  the  wall  normal  distance  of  up  to  1  mm.  This  comparison  shows  a 
boundary-layer  calculated  by  STABL  that  is  approximately  1  mm  thinner  than  that 
measured  by  Rufer.  The  maximum  difference  in  the  profile  shape  is  approximately 
20%  with  respect  to  the  massflux  scale.  These  measurements  are  being  repeated 
using  the  new  system  to  try  to  resolve  these  differences. 


88 


Figure  5.8:  Amplification  rate  curves  for  Rufer’s  blunt  cone  at  125  psia.  The  spacing 
between  the  curves  is  two  inches. 


s(m) 


Figure  5.9:  N  factor  curves  for  Rufer’s  blunt  cone  at  125  psia.  The  spacing  between 
the  curves  is  20  kHz. 
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co  (kHz) 

Figure  5.10:  Frequency  spectra  measured  by  Rufer  at  x  =  14.25  inches  on  the  blunt 
cone  at  125  psia. 


Figure  5.8  shows  amplification  rate  curves  for  every  two  inches  in  the  range 
4  in  <  x  <  18  in.  Well  defined  second  mode  peaks  are  present.  The  bandwidth  of 
the  unstable  region  is  approximately  150  kHz  near  the  nose  and  100  kHz  farther  aft. 
Figure  5.9  shows  N  factor  curves  for  this  case.  The  computation  was  carried  out 
with  10  kHz  spacing,  but  only  20  kHz  spacing  is  shown  in  the  figure  for  clarity.  A 
clear  envelope  is  evident.  N  =  5.5  is  first  reached  by  a  240  kHz  wave  at  s  =  0.388 
m  (x  =  15.2  in). 

Figure  5.10  shows  frequency  spectra  measured  by  Rufer  at  x  =  14.25  inches.  A 
sharp  peak  is  evident  near  the  boundary-layer  edge  with  no  corresponding  peak  near 
the  wall.  The  peak  frequency  of  approximately  240  kHz  is  close  to  the  frequency  with 
the  highest  N  factor  calculated  by  STABL,  which  is  260  kHz  at  that  location.  The 
frequency  calculated  by  STABL  to  have  the  highest  amplification  rate  at  x  =  14.25 
inches  is  226  kHz. 
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Figure  5.11:  A  comparison  of  the  mass  flux  calculated  by  STABL  with  the  uncali¬ 
brated  hot  wire  measurements  of  Rufer  at  two  locations  on  the  blunt  cone  at  45  psia. 


5.2.2  Blunt  Cone:  45  psia 

Figure  5.11  shows  a  comparison  of  the  boundary  layer  profiles  computed  by 
STABL  and  measured  by  Rufer  using  the  telescope  and  lens  system  for  the  blunt 
cone  at  45  psia.  No  discernable  difference  is  present  in  the  computed  and  measured 
boundary  layer  thicknesses.  In  addition,  the  shape  of  the  profiles  agrees  well,  with 
a  maximum  deviation  of  any  point  of  less  than  5%.  The  same  scale  is  used  for  both 
sets  of  experimental  data. 

Figure  5.12  shows  amplification  rate  curves  calculated  by  STABL  for  several 
axial  locations.  The  trend  here  is  slightly  different  than  that  seen  for  the  sharp  and 
blunt  cones  at  125  psia,  in  that  the  largest  amplification  rate  seen  is  not  the  most 
forward  curve.  The  trends  seen  in  the  amplification  rate  plots  for  the  other  cases 
are  probably  caused  by  the  computations  not  being  carried  out  far  enough  forward 
to  capture  behavior  just  aft  of  the  critical  point.  This  should  not  affect  the  N  factor 
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Figure  5.12:  Amplification  rate  curves  for  Rufer’s  blunt  cone  at  45  psia.  The  spacing 
between  the  curves  is  two  inches. 


calculations,  since  the  high  frequency  disturbances  that  begin  to  amplify  first  damp 
out  well  before  the  N  factor  envelope  reaches  the  transition  value.  The  N  factor 
calculations  were  begun  far  upstream  of  where  the  transition-causing  modes  first 
begin  to  amplify. 

Figure  5.13  shows  N  factors  for  the  blunt  cone  at  45  psia.  Figure  18  in  Rufer 
and  Schneider  [64]  shows  suspected  second-mode  instabilities  at  x  =  18  inches  and 
x  =  20  inches  but  not  at  x  =  22  inches,  suggesting  transition  occurred  between 
20  and  22  inches.  Using  the  curves  in  Figure  5.13,  these  locations  correspond  to 
N  factors  of  3.6  and  3.7.  Improved  frequency  resolution  in  the  computation  might 
increase  the  latter  number  to  approximately  4.0.  The  peaks  of  the  power  spectra  in 
Figure  18  of  [64]  occur  at  approximately  150  kHz  and  130  kHz.  The  most  amplified 
frequencies  in  Figure  5.13  at  those  locations  are  140  kHz  and  130  kHz,  respectively. 
Agreement  is  again  within  10%. 
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Figure  5.13:  N  factors  for  Rufer’s  blunt  cone  at  45  psia.  The  spacing  between  the 
curves  is  10  kHz. 


Appendix  D  provides  all  of  the  input  conditions  needed  to  obtain  the  mean  flow 
for  this  case  and  the  stability  results  shown  in  Figure  5.13.  Numerical  output  from 
selected  points  is  also  tabulated. 

5.3  Sensitivity  Studies 

In  order  to  gauge  the  effect  of  the  various  assumptions  made  in  setting  up  the 
boundary  conditions,  several  cases  were  analyzed  with  slight  perturbations  to  the 
boundary  conditions.  The  blunt  cone  at  45  psia  was  chosen  as  the  baseline  case 
due  to  the  available  high  quality  experimental  data,  the  fact  that  all  flight  vehicles 
have  some  degree  of  bluntness,  and  the  relative  ease  of  numerical  solution  of  blunt 
relative  to  sharp  cones.  The  test  conditions,  summarized  in  Table  5.3,  were  chosen 
to  represent  a  reasonable  band  of  uncertainty  in  the  nominal  conditions.  Each  case 
changed  only  one  variable  from  the  baseline  values. 
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(a)  Case  A  (b)  Case  B 


(c)  Case  C 


(d)  Case  D 


(e)  Case  E  (f)  Case  F 


Figure  5.14:  N  factor  curves  obtained  with  various  perturbations  to  nominal  bound¬ 
ary  conditions  for  Rufer’s  45  psia  blunt  cone. 
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Table  5.3:  Frequencies,  N  factors,  and  transition  locations  for  a  variety  of  input 
conditions  for  Rufer’s  blunt  cone  at  a  nominal  driver  tube  pressure  of  45  psia.  Fre¬ 
quencies  and  N  factors  are  for  x  =  20  inches. 


Case 

M oo 

Po  (psia) 

(K) 

Tw  (K) 

lo  (kHz) 

N 

Sn= 2.5  (mm) 

Base 

5.8 

45 

56 

300 

130 

3.63 

340 

A 

5.7 

45 

56 

300 

130 

3.47 

338 

B 

5.9 

45 

56 

300 

130 

3.47 

345 

C 

5.8 

40 

56 

300 

120 

3.19 

379 

D 

5.8 

35 

56 

300 

110 

2.75 

438 

E 

5.8 

45 

433 

300 

130 

2.82 

436 

F 

5.8 

45 

56 

distl 

130 

3.61 

339 

Figure  5.15:  N  factors  and  transition  locations  obtained  with  various  perturbations 
to  nominal  boundary  conditions  for  Rufer’s  45  psia  blunt  cone. 


95 


Figure  5.14  shows  the  N  factor  curves  for  the  six  perturbed  cases,  and  the  re¬ 
sults  are  summarized  in  Table  5.3.  Figure  5.15  graphically  depicts  the  N  factors 
and  transition  arc  lengths  shown  in  Table  5.3.  The  most  amplified  frequency  and 
maximum  N  factor  are  taken  at  x  =  20  inches  ( s  =  0.5121  m).  A  transition  arc 
length  is  also  given  based  on  an  Ntr  of  2.5.  Although  this  is  considerably  smaller 
than  the  expected  transition  N  factor,  this  was  the  highest  round  number  that  was 
reached  by  all  of  the  cases.  In  future  studies  of  this  type,  the  computational  domain 
should  be  extended  to  ensure  that  a  reasonable  transition  N  factor  is  reached  in  the 
domain.  Although  this  N  factor  is  low,  the  envelope  curve  is  reasonably  linear,  and 
the  trends  should  still  hold. 

Cases  A  and  B  represent  a  small  difference  in  the  freestream  Mach  number.  M <*,  is 
a  function  of  the  area  ratio  in  the  tunnel,  which  changes  with  p0.  Pitot  measurements 
with  the  cone  in  the  M6LT  showed  an  approximate  range  of  5.7  <  <  5.9,  but 

values  beyond  the  extremes  were  also  seen.  A  2%  change  in  M <*,  caused  a  4%  change 
in  the  N  factor  at  x  =  20  inches  and  a  1%  change  in  the  location  where  N  =  2.5. 
The  stability  results  do  not  appear  overly  sensitive  to  the  freestream  Mach  number. 

Cases  C  and  D  represent  a  change  in  the  tunnel  stagnation  pressure.  The  amount 
is  representative  of  the  typical  pressure  drop  during  the  course  of  the  run.  All  of 
Rufer’s  instability  measurements  were  taken  within  the  first  second  of  the  run  to 
minimize  the  effects  of  this  pressure  drop.  A  22%  decrease  in  p0  caused  a  24% 
decrease  in  N  and  a  29%  increase  in  s.  In  addition,  this  was  the  only  variable  which 
caused  a  significant  change  in  the  most  unstable  frequency.  Since  this  perturbation 
is  not  an  uncertainty,  but  rather  a  natural  change  in  the  test  conditions,  the  best 
results  will  be  obtained  if  all  measurements  are  taken  at  the  same  time,  as  is  currently 
done,  rather  than  using  the  full  test  time. 

Case  E  is  the  result  obtained  when  the  freestream  vibrational  temperature  is  set 
equal  to  the  stagnation  temperature.  As  discussed  earlier,  if  the  flow  in  a  wind  tun¬ 
nel  expands  quickly  enough,  the  translational  and  vibrational  energy  modes  may  not 
have  time  to  equilibrate.  The  stagnation  temperature  in  the  M6LT  is  almost  50% 
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Figure  5.16:  Wall  temperature  distribution  “distl”  used  in  sensitivity  study  on 
Rufer’s  45  psia  blunt  cone. 


lower  than  the  threshhold  value  proposed  by  Bertolotti  [57]  of  800  K.  In  addition, 
the  very  long  nozzle  in  the  M6LT  makes  flow  relaxation  more  likely  than  in  a  con¬ 
ventionally  designed  tunnel.  However,  some  degree  of  thermal  nonequilibrium  may 
be  present,  and  this  study  will  assess  the  potential  effect.  The  baseline  case  assumes 
full  equilibrium,  while  case  E  assumes  frozen  conditions.  The  effect  of  the  frozen 
conditions  is  a  22%  decrease  in  N  and  a  28%  increase  in  s.  This  is  a  significant 
effect,  and  a  computation  of  the  flow  of  the  tunnel  should  be  made  to  determine  the 
degree  of  thermal  nonequilibrium  present. 

Case  F  was  obtained  by  assuming  the  wall  temperature  had  the  distribution 
shown  in  Figure  5.16.  The  stagnation  point  wall  temperature  was  set  close  to  the 
freestream  stagnation  temperature,  and  the  wall  temperature  was  assumed  to  drop 
rapidly  and  fall  below  room  temperature  by  the  aft  end.  The  actual  distribution 
was  arbitrary,  but  this  was  assumed  to  be  a  worst-case  scenario  with  respect  to  the 
difference  from  a  constant  300K  wall.  The  effect  of  this  distribution  is  negligible, 
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Figure  5.17:  Amplification  rate  curves  at  x  =  20  inches  for  baseline  and  perturbed 
boundary  conditions. 


with  less  than  a  1%  effect  on  both  N  and  s.  This  suggests  that  additional  effort 
does  not  need  to  be  spent  to  determine  more  precisely  the  actual  wall  temperature 
distribution. 

The  amplification  rate  curves  at  x  =  20  inches  for  the  baseline  and  each  of  the 
perturbed  conditions  are  shown  in  Figure  5.17.  A  clear  trend  of  decreasing  frequency 
and  peak  amplification  rate  is  evident  with  increasing  Mach  number.  As  with  the  N 
factors,  the  largest  difference  from  the  baseline  case  is  seen  with  cases  C  and  D,  which 
correspond  to  the  total  pressure  drop.  This  has  the  effect  of  thickening  the  boundary 
layer,  which  decreases  the  unstable  frequency.  Interestingly,  although  freestream 
thermal  nonequilibrium  in  case  E  has  a  large  effect  on  the  N  factors  compared  to  the 
other  cases,  the  effect  on  the  local  amplification  rate  is  small.  Almost  no  difference 
is  seen  for  the  different  wall  temperature  in  case  F. 
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6.  Conclusions 

Future  systems  involving  hypersonic  flight  will  require  better  methods  of  laminar 
to  turbulent  boundary-layer  transition  prediction.  The  current  methods,  which  in 
most  cases  are  empirical  correlations  combined  with  extrapolation  from  wind-tunnel 
tests,  are  inadequate  to  reduce  design  margins  and  expand  the  state  of  the  art. 
The  semi-empirical  eN  method  coupled  with  LST  or  the  PSE  shows  promise  for 
some  conditions,  but  these  methods  have  not  yet  successfully  transitioned  from  the 
research  to  the  design  setting.  The  STABL  code  package  is  designed  to  make  that 
transition,  but  before  it  can  be  used  by  researchers  and  accepted  by  designers,  it  will 
require  extensive  verification  and  validation. 

STABL  and  its  stability  code,  PSE-Chem,  are  currently  transforming  from  an 
arcane  research  code  designed  for  a  single  type  of  problem  to  a  user-friendly,  broadly- 
applicable  prediction  tool.  As  the  first-ever  user  from  outside  the  developer’s  research 
group,  the  author  has  been  instrumental  in  suggesting  and  testing  improvements  to 
the  user  interface,  method  of  operation,  and  documentation.  Several  bugs  with 
varying  degrees  of  severity  have  been  uncovered.  The  experiences  of  the  author  have 
helped  STABL’s  developers  to  understand  the  areas  that  need  to  be  improved  to 
make  the  tool  accessible  enough  for  design  use.  Several  cases  have  been  analyzed  in 
an  effort  to  verify  and  validate  STABL. 

The  Stetson  et  al.  [33]  experiment  is  the  best-known  example  of  second  mode  wave 
growth  on  a  blunt  cone.  The  case  has  been  computed  by  many  other  researchers, 
which  makes  it  an  ideal  test  case  for  code  verification.  Many  issues  were  addressed 
to  improve  the  agreement  between  the  results  of  STABL  and  those  of  the  other 
researchers. 

This  effort  exposed  the  inadequacy  of  the  method  originally  implemented  in 
STABL  for  computing  the  viscosity  of  low-temperature  flows,  such  as  those  com- 
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monly  found  in  hypersonic  wind  tunnels.  Modifications  that  were  made  to  use  a 
blended  viscosity  law  improved  agreement  with  experimental  data  and  helped  to  ex¬ 
tend  the  applicability  of  STABL  to  the  low-temperature  regime.  A  related  issue  was 
the  use  of  the  proper  freestream  density.  Many  CFD  codes  require  a  unit  Reynolds 
number  input,  but  the  unit  Reynolds  number  reported  in  Stetson  et  al.  is  not  consis¬ 
tent  with  the  other  reported  quantities  if  a  Sutherland  viscosity  law  is  used.  Thus, 
any  researcher  that  used  the  unit  Reynolds  number  specified  by  Stetson  et  al.  as 
an  input  to  a  code  that  uses  the  Sutherland  law  modeled  the  flow  conditions  of  the 
experiment  incorrectly.  In  addition,  the  boundary-layer  edge  detection  algorithm 
was  modified  to  a  more  robust  scheme  which  allows  it  to  give  an  accurate  location 
even  when  overshoots  in  the  total  enthalpy  profile  exist.  The  sensitivity  of  several 
other  factors,  including  grid  convergence,  wall  temperature,  user-specified  options, 
freestream  thermal  non-equilibrium,  normalization,  compilers,  and  PSE-Chem  nu¬ 
merical  behavior  were  all  also  examined. 

The  transition  location  predicted  by  STABL  for  the  3.81  mm  Stetson  cone  agreed 
well  with  that  predicted  by  other  researchers.  The  location  where  N  =  5.5  is  within 
about  5%  of  the  other  values  in  the  literature.  However,  when  LST  is  used  to  cal¬ 
culate  the  amplification  rate  at  s/rn  =  175,  STABL  calculates  a  peak  amplification 
rate  that  is  more  unstable  than  any  of  the  data  sets  used  for  comparison.  The  differ¬ 
ence  is  approximately  equal  to  the  scatter  in  the  other  results.  Furthermore,  when 
the  viscosity  model  and  method  of  calculating  the  freestream  density  were  changed, 
the  results  calculated  by  STABL  shifted  by  amounts  similar  to  the  overall  scatter 
in  the  results.  This  shows  a  high  sensitivity  to  the  physical  models  employed,  and 
work  should  continue  in  validating  these  models  for  various  flow  regimes.  However, 
this  viscosity-model  uncertainty  is  currently  less  than  the  overall  uncertainty  in  the 
eN  method,  so  this  is  not  a  critical  issue  at  this  point. 

Computations  were  also  conducted  to  verify  the  ability  of  STABL  to  predict  first 
mode  disturbances.  The  maximum  N  factor  at  the  transition  location  was  12.2, 
which  is  slightly  higher  than  expected.  The  N  factor  computed  for  a  single  mode 
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differed  by  30%  compared  to  a  new  corresponding  calculation  using  the  eMahk  code. 
Analysis  of  the  eMahk  mean  flow  using  PSE-Chem  suggests  that  the  differences  could 
be  caused  by  differences  between  the  DPLR2D  Navier-Stokes  solution  and  the  eMahk 
similarity  solution. 

Analysis  of  an  adiabatic  large-bluntness  cone  at  Mach  8  showed  first-mode  am¬ 
plification  over  most  of  the  body.  Second-mode  amplification  began  farther  aft  at  an 
edge  Mach  number  of  5.7,  but  its  amplitude  rapidly  overtook  that  of  the  first  mode 
and  became  dominant  at  Me  =  6.0.  Agreement  with  the  results  of  Rosenboom  et 
al.  [38]  was  excellent  for  first  mode  transition  location  and  most  unstable  wave  an¬ 
gle,  although  there  was  a  significant  difference  in  frequency.  Second  mode  transition 
frequencies  matched  to  within  3%,  and  the  transition  location  agreed  to  within  10%. 

Computations  were  also  made  to  match  the  ongoing  hot-wire  experiments  of 
Rufer  in  the  Purdue  Mach  6  Ludwieg  Tube.  These  experiments  should  provide  an 
additional  source  of  badly-needed  validation  data.  Good  agreement  in  the  boundary 
layer  shape  and  thickness  is  obtained  for  both  blunt  and  sharp  cones  when  the  tele¬ 
scopic  height-measurement  system  is  used  in  the  experiments.  Agreement  between 
experiment  and  computations  in  the  most  unstable  frequency  is  within  10%  for  all 
cases  shown  in  Rufer  and  Schneider  [64],  Transition  N  factors  based  on  locations 
inferred  from  hot-wire  spectra  ranged  from  3.6  to  greater  than  5.4.  This  matches  the 
expected  trend  that  the  transition  N  factors  will  be  lower  than  in  Tunnel  B,  since 
the  smaller  test  section  means  that  the  noise  intensity  will  be  higher.  The  transition 
N  factors  decrease  as  the  tunnel  stagnation  pressure  decreases.  This  matches  the 
unit  Reynolds  number  effect  observed  in  conventional  wind  tunnels,  which  says  that 
transition  Reynolds  number  decreases  as  unit  Reynolds  number  decreases  [10]. 

Sensitivity  studies  were  also  conducted  to  gauge  the  effect  of  various  assump¬ 
tions  that  were  made.  A  dependence  on  freestream  Mach  number  was  observed,  as 
a  change  of  2%  in  the  Mach  number  caused  a  4%  change  in  the  N  factor  at  a  given 
location.  The  normal  change  in  the  driver  tube  pressure  over  the  course  of  a  tunnel 
run  caused  a  25%  change  in  the  N  factor  and  a  15%  change  in  the  frequency.  This 
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reiterates  the  importance  of  taking  measurements  during  the  same  driver-tube  pres¬ 
sure  window  for  each  run.  Freestream  thermal  non-equilibrium  was  found  to  have  a 
potentially  large  effect  on  the  stability  computations,  with  a  frozen  freestream  vibra¬ 
tional  temperature  resulting  in  a  22%  decrease  in  the  N  factor,  although  the  actual 
degree  of  freestream  thermal  nonequilibrium  is  unknown.  Almost  no  effect  was  seen 
from  a  different  wall  temperature  distribution. 

A  general  observation  is  better  agreement  is  found  in  overall  transition  locations 
computed  using  different  computational  methods  than  in  the  individual  computa¬ 
tions  used  to  create  the  overall  prediction.  This  was  seen  in  Stetson’s  3.81  mm  blunt 
cone,  where  STABL  calculated  a  local  amplification  rate  that  was  higher  than  other 
researchers’  results,  but  the  predicted  transition  location  fell  within  the  scatter.  It 
was  also  seen  for  the  sharp  cone  at  Mach  3.5.  STABL  and  eMahk  calculated  N  fac¬ 
tors  for  a  single  mode  that  disagreed  by  30%,  but  the  transition  N  factor  STABL 
calculated  was  within  15%  of  the  expected  range  for  low  disturbance  environments. 
Since  the  N  factor  is  an  integrated  quantity,  and  determining  the  transition  location 
through  the  N  factor  envelope  involves  using  many  individual  growth  rates,  it  seems 
logical  that  this  process  would  tend  to  smooth  out  differences  in  the  growth  rates. 

It  is  interesting  that  the  N  factors  seen  for  the  sharp  cone  at  Mach  3.5  were  larger 
than  expected,  whereas  those  seen  for  Rufer’s  cases  were  smaller  than  expected.  It 
is  important  to  remember  that  the  eN  method  is  semi-empirical.  The  N  factor  that 
best  predicts  transition  is  dependent  on  freestream  noise  and  other  factors,  and  it 
needs  to  be  chosen  from  comparisons  to  experimental  data.  When  Malik  [61]  first 
computed  the  sharp  cone  case,  he  appears  to  have  performed  a  total  of  seven  N  factor 
calculations.  In  contrast,  twenty-one  years  of  advances  in  hardware  and  software 
allowed  calculation  of  480  individual  N  factors  here.  It  does  not  seem  unreasonable 
that  this  could  result  in  a  different  transition  N  factor. 

A  great  deal  of  progress  has  been  made  in  transforming  STABL  from  a  research 
code  to  a  design  tool.  The  computations  presented  here  help  to  provide  confidence 
in  the  ability  of  STABL  to  analyze  boundary-layer  stability  and  predict  transition. 
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These  results  can  also  serve  as  a  benchmark  to  allow  other  researchers  using  STABL 
to  verify  their  method  of  operation.  Additional  computations  on  different  cases  in  a 
variety  of  flow  regimes  should  be  performed  to  further  this  effort. 
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Appendix  A:  Typical  N-Factor  Script 

An  example  script  used  to  run  multiple  combinations  of  starting  location,  fre¬ 
quency,  and  spanwise  wavenumber  is  given  below. 

# ! /usr/bin/perl 

#  Reference  the  pse_scripting  library  routines 
require  . /lib/pse_scripting . lib" ; 

#  Set  the  directory  in  which  to  store  the  results 
my  $od=  "/home/trobarge/pseruns/133" ; 

my  $pid=  $$; 

#  Set  number  of  processors 
my  $np=2; 

#  Set  up  the  parameter  space  for  the  test  matrix 
my  @sdists=  map  {$_/50;}  ( 1 . .3) ; 

my  @freqs=  map  {$_*25;}  (1..10); 
my  @betas=  map  {$_*1000;}  (1..8); 

#  Set  PSE-Chem  options  to  be  used  for  all  runs 

my  °/0params=  (iest_omeg=>0 ,  iftest3=>l,  lstest3=>0, 
flow=>0.2,  fhigh=>1.0,  istep_march=>l , 
mean_f low_f ile=>" J . . /out/f irstl45 .bin’ " , 
ichem_on=>l,  ivib_on=>l,  idiff_on=>l) ; 

#  Create  status  file  -  deletion  stops  the  run 

my  $statusfile=  psechem_statusf ile_new(" ’ continue ’ ) ; 
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print  " —  Status  file  is:  $statusf ile\n" ; 

#  Make  the  main  output  directory  if  it  doesn't  already  exist 
if  ( ! -d  $od)  { 

mkdir($od) ; 
if  ( ! -d  $od)  { 

die  "Unable  to  make  output  directory  ’ $od'\n"; 

> 

> 

#  Loop  over  frequencies  and  set  options  for  each  loop 
foreach  my  $f  (@freqs)  { 

$params{omegal}=$f ; 

$params{omega2}=$f ; 

#  Loop  over  starting  locations  and  set  option  for  each  loop 
foreach  my  $s  (@sdists)  { 

$params{"start_sdist"}=  $s; 

#  Loop  over  betas  and  set  options  for  each  loop 

foreach  my  $beta  (©betas)  { 

$params{betaO_local}=$beta ; 

$params{betaO_pse}=$beta; 

#  Name  run  according  to  individual  parameters 

my  $run="$f-$s-$beta" ; 

#  Options  to  pass  to  PSE-Chem  solver 

my  @oargs=  (-np=>$np,  -copy_to=>$od,  -id=>"$run". 


-args=>">  $od/psechem.$run.out" , 
-outdir=>" . . /out/pserun . $pid . data") ; 
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#  Pass  options  to  input  file  and  solver  and  run  using  library  routine 

print  " —  Running  PSE-Chem  on  case  $run\n" ; 

my  %ret=  psechem_run_solver(-inpf  ile=>Y/0params , @oargs) ; 

psechem_check_error(°/0ret)  ; 

#  Check  status  file  to  see  whether  to  continue 

my  $status=  psechem_statusf ile_check($statusf ile) ; 
if  ($status=~/abort/)  { 

print  "Stopping.  Status  file  contains:  $status\\\n" ; 
if  (-f  $statusfile)  {unlink ($statusf ile) ;} 
exit ; 

> 

> 

> 

} 

#  Remove  status  file  and  exit  when  all  runs  finished 
if  (-f  $statusfile)  {unlink ($statusf ile) ;} 

print  "\n —  Done  — \n\n" ; 
exit ; 
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Appendix  B:  STABL  User’s  Guide 

This  appendix  is  designed  to  be  a  guide  to  using  STABL  from  a  user’s  perspective. 
It  supplements  the  extensive  documentation  included  with  the  STABL  distribution. 
This  is  believed  to  be  accurate  as  of  July  2005,  but  since  STABL  is  still  under 
active  development,  many  of  the  details  provided  here  may  become  obsolete,  and 
the  STABL  documentation  should  be  considered  the  authoritative  source. 

The  opinions  given  in  this  appendix  are  my  own,  and  they  should  not  be  construed 
as  those  of  STABL’s  developers.  I  will  describe  the  lessons  I  have  learned,  the  way 
I  used  STABL  and  set  up  the  environment,  and  in  some  cases  the  way  I  wish  I  had. 
Individual  preferences  will  vary,  and  I  don’t  claim  to  have  found  the  best  way  to  do 
everything,  so  use  this  as  a  guide,  not  as  a  directive. 

B.l  Prerequisites 
Hardware 

STABL  parallelizes  very  efficiently,  allowing  it  to  effectively  utilize  the  resources 
of  a  large  cluster  supercomputer,  but  it  can  also  run  fully  on  a  single  workstation. 
This  section  will  summarize  my  experiences  and  recommendations. 

I  was  given  access  to  an  812  processor  cluster  at  Sandia  National  Laboratories, 
and  I  used  that  computer  to  run  most  of  my  mean  flows.  I  typically  used  8-16  pro¬ 
cessors  at  a  time,  and  grids  with  150,000  points  would  converge  in  approximately  six 
hours  with  no  user  intervention.  If  necessary,  I  could  be  more  aggressive  in  increas¬ 
ing  the  CFL  number  and  freezing,  and  significantly  reduce  the  time  to  converge,  but 
this  was  normally  not  worth  the  user-time  involved,  given  the  complexity  caused  by 
the  batch  system  and  the  need  for  remote  operation  and  file  transfers.  I  also  occa¬ 
sionally  used  the  Sandia  cluster  for  PSE-Chem  runs,  particularly  when  I  ran  scripts 
involving  large  numbers  of  global  and  local  frequencies,  since  those  portions  of  the 
code  parallelize  well. 
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I  used  a  dual  processor  Symmetric  Multi-Processing  (SMP)  workstation  for  the 
majority  of  my  time  with  STABL.  It  was  a  Dell  Precision  650,  which  was  their 
largest  workstation  model  at  the  time  it  was  purchased.  The  workstation  contained 
two  3.06  GHz  Intel  Xeon  processors,  each  with  1  MB  of  L3  cache.  The  system  had 
2  GB  of  ECC  SDRAM  and  an  80  GB  IDE  hard  drive.  Graphics  were  rendered  by 
a  dual-monitor-capable  nVidia  QuadroFX  500  graphics  card  with  128  MB  of  video 
RAM.  A  CD-RW/DVD  drive,  3.5  inch  floppy  drive,  keyboard,  and  mouse  were  also 
included.  No  monitors  were  purchased,  as  suitable  used  monitors  were  available  free 
of  charge  from  the  University.  The  total  cost  of  the  system  in  May  2004,  including 
academic  and  volume  pricing,  was  $2,628.80. 

In  retrospect,  all  of  the  computations  presented  in  this  thesis  could  have  been 
done  on  the  workstation  alone,  although  the  Rufer  sensitivity  studies  would  have 
taken  considerably  longer.  The  mean  flows  could  be  obtained  in  2-3  days,  and 
the  workstation  could  still  be  used  for  other,  less  processor-intensive  tasks  such  as 
data  analysis  in  the  mean  time.  There  would  be  a  significant  effect  on  the  system 
response  time  for  other  tasks,  but  if  the  full  capability  of  4  GB  of  RAM  had  been 
purchased  and  the  process  prioritizing  ability  of  Linux  was  used,  the  effect  could 
be  greatly  decreased.  However,  the  nature  of  the  work  in  this  thesis  was  limited 
by  analysis  time,  rather  than  CPU  cycles.  If,  for  example,  a  parametric  study  was 
performed  to  obtain  transition  Reynolds  numbers  for  many  combinations  of  flow  and 
geometry  variables,  the  analysis  load  for  each  individual  solution  would  be  much 
smaller,  and  CPU  cycles  would  become  the  limiting  factor.  Therefore,  the  necessity 
of  supercomputer  access  will  probably  depend  on  the  manner  of  usage  of  STABL. 

Irregardless  of  supercomputer  availability,  a  quality  workstation  is  essential  for 
intensive  use  of  STABL.  The  workstation  will  at  a  minimum  be  used  for  compilation, 
post-processing,  and  data  analysis.  I  found  the  workstation  described  above  to 
be  excellent  for  all  of  those  tasks.  A  high  end  graphics  card  is  not  necessary  for 
2D/axisymmetric  data  visualization,  and  the  reasonably  priced  card  chosen  allowed 
me  to  use  two  monitors  to  double  the  effective  size  of  my  workspace,  which  greatly 
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increased  my  productivity.  I  found  the  processor  speed  and  amount  of  RAM  fully 
adequate  for  my  purposes.  Changes  I  would  make  now  would  be  to  purchase  a  DVD 
writer  to  assist  in  data  backup  and  a  SCSI  hard  drive  instead  of  IDE  for  faster  data 
transfer. 

I  found  it  useful  to  have  a  fixed  IP  address  and  registered  host  name.  Since  I 
accessed  the  Sandia  cluster  remotely  through  SSH,  this  allowed  me  to  connect  back 
to  my  workstation  from  Sandia  to  transfer  files  in  a  single  step  rather  than  through 
a  proxy.  In  addition,  I  could  connect  to  the  workstation  from  home  in  the  evenings 
and  weekends  to  monitor  the  status  of  ongoing  computations. 

Software 

STABL  was  written  in  and  for  a  UNIX  environment.  Many  varieties  of  UNIX  are 
available;  I  recommend  using  Linux  if  a  new  system  is  used.  All  of  my  computations 
were  done  using  Red  Hat  Linux  9.  The  workstation  came  with  Windows  loaded; 
Linux  was  installed  locally,  with  Windows  retained  on  a  second  partition.  Other 
varieties  of  UNIX  should  work  fine  once  initial  setup  is  complete,  although  there 
may  be  some  extra  work  if  this  is  the  first  time  STABL  has  been  used  on  that 
particular  variety.  It  may  be  possible  to  run  STABL  under  Windows  using  Cygwin, 
which  is  a  freely  available  program  that  provides  a  UNIX  environment  on  a  Windows 
operating  system.  Currently,  this  has  never  been  tried  extensively,  and  I  would  highly 
recommend  Windows  users  either  use  a  second  computer  or  create  a  dual-boot  system 
with  Linux. 

STABL  is  distributed  in  source  form,  so  C  and  Fortran  compilers  are  needed  to 
create  the  binary  executable  files.  I  originally  used  the  GNU  compilers,  but  later 
switched  to  the  Portland  Group  compilers  for  reasons  described  in  Section  3.2.9.  The¬ 
oretically,  the  Intel  compilers  should  work  also,  although  these  have  not  been  used 
yet  to  my  knowledge.  The  Portland  Group  C  and  Fortran  compilers  were  down¬ 
loaded  from  the  website  and  installed  on  the  workstation  using  the  license  manager 
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option.  This  was  slightly  more  complicated  to  set  up,  but  it  allowed  multiple  users 
access  to  the  compilers.  The  details  are  described  in  the  compiler  documentation.  A 
Perl  interpreter  is  also  required  for  the  GUI,  scripting  capabilities,  and  many  of  the 
additional  utilities.  This  is  usually  installed  by  default  on  UNIX  systems,  and  a  Win¬ 
dows  version  can  be  downloaded  free  of  charge  from  http://www.activestate.com. 
The  STABL  documentation  contains  detailed  information  on  the  installation  of  ad¬ 
ditional  Perl  modules  that  are  required. 

DPLR2D  and  PSE-Chem  use  MPI  libraries  for  message  passing  between  multi¬ 
ple  processors.  Several  versions  of  MPI  are  available;  I  installed  the  MPICH  MPI 
libraries,  which  are  freely  available  online.  The  LAM  MPI  libraries  were  also  tested, 
and  no  effect  was  seen  on  the  results.  I  found  it  best  to  set  up  the  compilers  be¬ 
fore  setting  up  the  MPI  libraries  so  that  the  mpif  90  and  mpicc  commands  could  be 
properly  configured.  The  Portland  Group  website  FAQ  contains  information  about 
setting  up  their  compilers  with  MPICH. 

The  commercial  data  visualization  software  Tecplot  was  used  extensively  for  anal¬ 
ysis  of  the  results.  STABL  is  designed  to  format  its  data  files  to  be  Tecplot-ready, 
although  it  would  not  be  very  difficult  to  change  the  source  code  so  that  another 
formatting  is  used.  A  single  commercial  license  of  Tecplot  currently  costs  $1600, 
and  an  educational  license  is  approximately  $800.  Many  institutions  will  have  a  site 
license  which  would  make  Tecplot  available  at  no  additional  cost,  as  was  the  case  for 
me. 

Knowledge 

STABL  is  an  advanced  system,  and  users  more  experienced  with  CFD  and  tran¬ 
sition  prediction  will  have  an  easier  learning  curve  than  others.  In  particular,  a 
working  knowledge  of  CFD  principles,  as  can  be  gained  in  an  introductory  course, 
will  be  very  helpful.  A  user  should  be  reasonably  proficient  in  using  the  command 
line  in  UNIX.  Knowledge  of  Fortran  is  only  strictly  necessary  to  understand  the 
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source  code,  but  that  would  be  beneficial,  since  the  code  is  still  under  active  devel¬ 
opment.  A  basic  knowledge  of  scripting  in  Perl  is  necessary  to  utilize  the  advanced 
operation  methodologies  currently  being  developed.  While  only  basic  knowledge  in 
these  areas  is  necessary  to  use  STABL,  more  advanced  users  will  be  better  able  to 
exploit  the  full  capabilities  of  the  software. 

B.2  Method  of  Operation 

The  STABL  documentation  gives  a  detailed  tutorial  to  allow  a  user  to  run  a  case 
with  all  of  the  correct  options  in  place.  While  this  is  sufficient  to  teach  the  new  user 
how  to  execute  the  proper  commands,  it  does  not  teach  why  those  options  were  set 
and  when  a  user  will  want  to  change  the  defaults.  This  section  is  an  attempt  to  fill 
that  gap  with  some  of  the  lessons  I  learned  over  the  course  of  my  time  using  STABL. 

Case  Definition 

The  first  step  in  analyzing  is  a  new  case  is  to  gather  the  relevant  information 
about  that  case.  What  you  need  depends  on  your  goals  for  the  analysis,  but  at  a 
minimum  you  will  need  the  model  geometry,  the  freestream  properties,  and  the  wall 
temperature  boundary  condition.  The  information  needed  for  the  model  geometry 
depends  on  the  class  of  problem.  For  a  blunt  cone,  you  will  need  the  nose  radius, 
the  half  angle,  and  the  arclength  from  the  sphere-cone  tangency  point  to  the  end. 
For  the  input  conditions,  you  will  need  the  freestream  Mach  number  or  velocity,  the 
freestream  static  translational  temperature  T^,  the  freestream  static  vibrational 
temperature  Tv,*,,  and  the  freestream  density.  You  will  most  likely  need  to  calculate 
the  density  from  other  parameters.  Be  careful  in  doing  this  to  ensure  you  are  doing 
it  in  a  way  that  is  consistent  with  the  goal  of  the  analysis;  see  Section  3.2.2  for  more 
discussion  of  this.  For  the  wall  temperature  condition  you  can  specify  a  constant 
temperature,  a  temperature  distribution,  radiative  equilibrium,  or  an  adiabatic  wall. 
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You  may  need  other  information  to  make  use  of  the  STABL  results.  If  you  want 
to  compare  to  experiments,  you  will  need  data  such  as  the  experimental  transition 
location  or  the  wall  temperature  or  pressure  distributions.  For  transition  prediction, 
you  will  need  to  specify  an  N  factor  appropriate  for  the  environment.  It  would  be 
helpful  to  know  the  relative  uncertainty  in  the  various  input  quantities  so  that  at 
least  basic  sensitivity  studies  can  be  performed. 

Grid  Generation 

The  first  step  in  the  analysis  is  to  generate  a  grid.  I  used  the  STABL  GUI  to 
run  the  grid  generator.  The  model  geometry  is  used  to  define  the  inner  boundary. 
The  outer  boundary  should  be  tailored  to  match  the  angle  of  the  shock  as  closely 
as  possible.  At  least  five  or  six  grid  points  should  stand  between  the  shock  and  the 
domain  outer  boundary  to  ensure  the  shock  is  well  captured. 

The  problem  with  matching  the  grid  to  the  shock  shape  is  that  the  shock  shape 
is  not  known  at  the  outset.  To  get  around  this,  I  recommend  using  two  grids.  To 
make  the  first,  use  the  graphs  in  NACA  Report  1135  [65]  or  Billig’s  correlation  [66] 
to  obtain  an  approximate  outer  shock  angle.  Apply  a  generous  margin  of  error 
and  create  a  large,  coarse  grid  that  you  are  sure  will  completely  capture  the  flow. 
Run  the  solution  for  this  case,  which  due  to  the  small  problem  size  and  low  aspect 
ratios  should  converge  very  quickly.  Postprocess  the  solution  and  plot  Mach  number 
contours  in  Tecplot.  I  like  to  make  the  shock  the  only  contour  level  shown  to  enhance 
the  contrast.  Then  make  a  second  grid  that  will  better  match  the  shock.  Read  the 
grid  file  into  Tecplot  and  overlay  it  on  the  contours.  I  like  to  activate  only  the 
contour  layer  in  the  original  solution  and  only  the  mesh  or  boundary  in  the  new 
grid.  Now  go  back  to  the  GUI  and  make  another,  better  grid  that  uses  the  lessons 
learned.  Continue  doing  this  until  your  new  grid  matches  the  shock  shape.  I  typically 
use  around  ten  line  segments  in  constructing  the  outer  grid  shape  for  a  blunt  cone. 
Tecplot  macros  can  make  the  iterations  go  much  faster.  Increase  the  number  of 
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points  in  the  second  grid.  The  spacing  factors  may  have  to  be  changed  to  place 
sufficient  points  within  the  boundary  layer.  The  boundary-layer  thickness  can  be 
approximately  determined  from  the  first  solution,  as  well. 

After  the  second  grid  is  created,  save  your  work  and  run  the  mean  flow  solver. 
Be  sure  to  check  that  your  grid  still  matches  the  shock  well  and  the  spacing  within 
the  boundary  layer  is  sufficient  in  the  final  solution. 

Mean  Flow  Solution 

The  next  step  is  to  obtain  the  laminar  mean  flow.  The  STABL  documentation 
gives  a  good  description  of  how  to  run  the  solver.  I  will  focus  on  how  to  set  up  the 
problem  and  various  choices  that  need  to  be  made.  Most  of  the  work  is  in  setting 
up  the  input  file.  Brief  descriptions  of  each  of  the  options  are  given  at  the  bottom  of 
the  sample  file.  I  will  focus  on  the  ones  that  you  are  most  likely  to  need  to  change 
from  the  defaults. 

The  grid  filename  is  self-explanatory,  iwall  should  be  1  for  an  adiabatic  wall,  2 
for  an  isothermal  wall,  3  for  a  specified  distribution,  or  4  for  radiative  equilibrium. 
For  an  isothermal  wall,  the  temperature  is  specified  by  Twall,  and  for  a  specified 
distribution  it  is  given  in  the  file  declared  as  the  wall  temperature  filename, 
istop  gives  the  total  number  of  iterations  to  perform  in  the  run.  If  you  restart  from 
a  previous  run,  the  iteration  count  resets  for  the  purposes  of  istop.  To  restart, 
change  iconr  to  1,  ensure  that  the  <BASE>.flow  file  is  in  the  .  ./out  directory,  and 
run  the  solver  with  the  normal  command. 

The  kbl  parameter  is  most  likely  to  be  overlooked.  In  layman’s  terms,  this  tells 
the  solver  at  which  j  value  it  should  toggle  a  limiter  that  results  in  an  increase  in  the 
effective  viscosity.  This  should  be  set  to  a  j  value  such  that  is  at  every  i  beyond  the 
boundary  layer  edge  but  inside  of  the  shock.  A  suitable  value  should  be  determined 
when  the  second  grid  is  created.  If  this  is  set  to  a  value  within  the  boundary  layer, 
discontinuities  may  appear  in  the  profile  derivatives. 
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nplot  specifies  the  interval  of  iterations  between  writing  restart  files.  I  adjust 
this  based  on  the  problem  size  and  system  used  so  that  it  writes  a  restart  file  every 
minute  or  two.  igeom  should  be  set  to  1  for  a  blunt  cone  or  3  for  a  sharp  cone,  ichem 
and  ivib  settings  are  at  your  discretion.  They  should  probably  both  be  1  unless 
you  are  trying  to  compare  to  another  computation  or  assess  the  real  gas  effects. 
The  initial  CFL  number  is  set  through  cf  1,  and  subsequent  numbers  to  change  to 
are  listed  after  the  mass  fractions.  The  CFL  is  changed  every  ijump  iterations. 
The  freestream  conditions  and  the  mass  fractions  will  also  need  to  be  set;  those  are 
self-explanatory. 

I  prefer  to  run  new  cases  on  my  workstation  first,  where  I  can  easily  monitor  the 
progress,  start  and  stop  easily,  and  experiment  with  different  CFL  numbers.  Once  I 
have  gained  confidence  that  the  solution  is  converging  well,  I  will  transfer  it  to  the 
cluster,  where  it  runs  faster  but  is  much  more  difficult  to  monitor. 

The  mean  flow  code  gives  you  the  option  to  change  the  freezing  level  and  CFL 
number  while  it  is  running.  This  capability  is  very  useful  in  a  batch-scheduling 
system  where  restarting  the  job  may  mean  going  to  the  end  of  the  queue.  However, 
if  you  increase  the  CFL  number  too  aggressively,  the  solution  may  blow  up.  For  this 
reason,  I  was  often  overly  conservative  and  would  let  the  solution  converge  slower 
than  it  needed  to.  However,  I  gained  the  flexibility  to  let  it  run  unattended  and  it 
would  eventually  converge. 

After  the  residual  has  leveled  off,  you  should  postprocess  it  and  look  for  several 
factors.  First,  plot  contours  of  the  residual  on  a  logrithmic  scale.  I  used  a  Tecplot 
macro  to  easily  automate  this.  I  always  saw  slightly  higher  residuals  near  the  nosetip 
in  the  boundary  layer,  but  other  “hot  spots”  should  be  examined  carefully.  Look  at 
the  shock  shape  to  make  sure  the  grid  matches  the  angle  while  maintaining  a  reason¬ 
able  standoff  distance.  A  good  check  is  to  probe  at  various  points  outside  the  shock 
to  ensure  that  the  freestream  values  were  maintained.  I  also  like  to  look  at  velocity 
and  temperature  profiles  to  ensure  that  they  look  reasonable  before  proceeding  to 
the  stability  analysis. 
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Stability  Analysis 

The  final  step  is  to  complete  the  stability  analysis.  This  step  is  particularly 
dependent  upon  the  goals  for  the  analysis.  Predicting  transition  is  relatively  easy. 
The  walkthrough  included  with  the  STABL  release  gives  good  examples  for  obtaining 
N  factors  and  LST  amplification  curves.  Almost  all  of  the  plots  I  show  in  this  thesis 
were  obtained  using  the  PSE  scripting  routines  with  setup  files  modified  from  the 
examples  provided.  To  obtain  N  factors  for  transition  prediction,  I  recommend 
running  a  two-parameter  test  matrix  with  combinations  of  frequency  and  starting 
location.  For  first  mode  N  factors,  the  j3  value  is  a  third  parameter  that  must  be 
varied. 

After  the  runs  are  completed,  use  command-line  wild-card  expansion  to  load 
all  of  the  psealpha  files  into  Tecplot.  Plot  each  of  the  N  factors,  and  select  the 
maximum  N  factors  for  each  frequency.  A  routine  has  just  been  added  to  STABL  to 
extract  the  most  unstable  modes  automatically,  but  I  have  not  tried  it  myself.  The 
predicted  transition  location  will  be  where  the  maximum  N  factor  first  reaches  the 
preset  value.  Amplification  rate  curves  and  eigenfunctions  can  be  plotted  if  desired. 
Other  work  will  be  determined  by  the  goals  of  the  study. 

There  are  several  input  file  parameters  that  I  change  on  a  regular  basis.  The  pa¬ 
rameters  mean_f low_f ile,  start_sdist,  iselect_freq,  betaO_local,  betaO_pse, 
and  iest_omeg  are  all  well  explained  in  the  documentation  and  in  the  descriptions  in 
the  input  file.  I  vary  num_coarse  and  num_f  ine  depending  on  the  application,  but 
I  rarely  increase  them  over  the  default  values.  I  normally  leave  ichem_on,  ivib_on, 
and  idiff_on  set  to  1  unless  I  am  specifically  testing  their  effect.  I  have  not  seen 
an  effect  between  the  various  settings  of  istr,  with  the  exception  of  the  uniform 
setting.  I  normally  leave  that  at  1,  since  I  use  high-quality  mean  flow  grids.  I  use  a 
bl_method  of  3,  but  you  should  check  to  see  if  the  boundary  layer  edge  location  is 
reasonable.  I  have  experimented  with  most  of  the  other  options,  and  I  recommend 
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that  other  users  do  the  same,  but  those  are  the  ones  I  generally  considered  when 
setting  up  a  new  case. 

B.3  Final  Thoughts 

Hopefully  this  users’  guide  has  filled  some  of  the  gaps  in  the  STABL  documen¬ 
tation  and  provided  a  slightly  different  perspective  on  the  working  of  the  code.  I 
would  encourage  you  to  ask  questions  and  to  try  out  the  many  options  available. 
The  PSE-Chem  source  code  is  well  documented  and  laid  out,  so  you  can  often  look 
directly  at  that  section  of  the  code  to  learn  the  differences  between  the  various  op¬ 
tions.  Take  careful  notes,  and  always  save  your  input  files.  If  you  have  your  input 
files,  you  can  recreate  your  data;  if  you  only  have  your  data,  you  may  not  be  able  to 
determine  how  you  obtained  it.  Best  wishes. 
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Appendix  C:  Boundary  Layer  Edge  Detection 

The  edge  detection  algorithm  was  modified  to  account  for  overshoots  in  the  total 
enthalpy  profiles.  The  file  STABL/cfd_solver/src/cfd_detect_bl_edge  .fpp  was 
the  only  file  modified.  A  partial  listing  of  the  source  code  is  shown  below. 

c  ***  Total  enthalpy  method 
if  (imethod . eq. 1)  then 

hp=  rbl_fact 

hf=  h(nj)  !  Freestream  value 

hw=  h(l)  !  Wall  value 

hd=  abs(hf  -  hw)  !  Difference 

c  ***  Initialize  variables 

jedge  =  0 
ioversht=  0 

c  ***  First  use  original  method 

do  j=  1 ,nj-l 

if  (abs(h(j)-hw)  .ge.  hp*hd)  then 

jedge  =  j 
jf irst=  j 
exit 
endif 
enddo 

if  (jedge. eq.O)  goto  999 

c  ***  Find  out  if  an  overshoot  exists 


122 


do  j  = 


c 

c 


enddo 


jf irst+1 ,nj 

if  ((h(j)-hw)  .gt.  hd/hp)  then  ! overshoot  exists 
ioversht=  1 
endif 

***  Next  line  makes  sure  we’re  past  overshoot  and  checks  the  flag 
***  to  make  sure  it  isn’t  fooled  by  random  fluctuations  in  the  freestream 
if  (iand(ioversht==l,h(j)  .It.  h(j-l)))  then 
if  ((h(j)-hw)  .It.  hd/hp)  then 
jedge=  j 
exit 
endif 
endif 
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Appendix  D:  Input  and  Output  Data  for  Rufer’s  45  psia  Blunt  Cone 

This  appendix  is  designed  to  provide  all  of  the  input  data  needed  to  recreate  the 
plot  shown  in  Figure  5.13.  The  case  is  Rufer’s  blunt  cone  with  a  stagnation  pressure 
of  45  psia.  Computations  were  performed  in  May  2005  using  STABL  1.29.2. 

The  mean-flow  input  file  is  provided  in  Section  D.l.  All  of  the  inputs  needed  for 
the  GUI  elliptical  nose  grid  generator  are  provided  in  Table  D.l.  A  script  that  can 
be  used  to  generate  the  PSE  data  is  printed  in  Section  D.2.  Mean-flow  output  data 
is  tabulated  for  several  locations  in  Table  D.2.  N  factors  for  several  frequencies  and 
locations  are  printed  in  Table  D.3. 

D.l  Mean  Flow  Input  File 

!  Grid  filename, 

’ . . /grids/grid_shannl50 . dat ’ 

!  chemistry  input  filename, 

’ . . /props/air_5sp_88.inpJ 
!  wall  temperature  filename, 

’ . . / inp/twall . inp ’ 

!  NASA  thermodynamic  properties  filename 
’ . . /props/nasa_1996_thermo .dat ’ 


itvd 

iorder 

i ext st 

kmax 

ivis 

iwall 

i2j 

-1, 

2, 

-1, 

4, 

1, 

2, 

2 

istop 

iplot 

iconr 

iaxi 

inor 

isn 

kbl 

50000, 

1, 

-1, 

1, 

U 

0, 

220 

nplot 

igrid 

ilt 

i2n 

igeom 

irm 

iej 

100, 

1, 

-1, 

5, 

1, 

1, 

-1 

ichem 

ivib 

itv 

itl 

irk 

iset 

icf ljmp 

1, 

1, 

1, 

-1, 

1, 

1 

20 
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machin 

density 

Tin 

Tvin 

Twall 

vin 

rvr 

5 . 80d0 , 

1 . 5046d-2 , 

5 . 6000dl , 

5 . 6000dl , 

3 . 00d2 , 

3. 06000d3, 

1 . 3d0 

cfl 

epsi 

epsj 

epsk 

pmul 

alpha 

yaw 

1.0d-3, 

0.3, 

0.3, 

0.3, 

0.5, 

O.OdO, 

O.OdO 

stime 

rconv 

rf reeze 

raccel 

vaccel 

emisw 

radrx 

10.0 

0.5d-7 

-1.0 

1 .  Od+5 

1 . Od+5 

0 . 8d0 

0 . 5d0 

!  Mass  fractions,  species  order:  N2,02,N0,N,0 
0.767000 
0.233000 
0.000000 
0.000000 
0.000000 

!  List  of  additional  CFL  numbers 
0.01 
0.1 


10. 

20. 

40. 

100.  ,5 
200.  ,5 
400.  ,  20 
1000. 


-1 
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Table  D.l:  Grid  input  conditions  for  Rufer’s  45  psia  blunt  cone. 


Major  (x)  nose  radius 

0.000508  m 

Minor  (y)  nose  radius 

0.000508  m 

Number  of  points  on  nose 

80 

Major  (x)  axis  upper  grid  radius 

.00141 

Minor  (y)  axis  upper  grid  radius 

.00120 

Upper  grid  x  offset 

0.0008 

Number  of  body-normal  points 

351 

Nose  grid  point  spacing  factor 

0.07 

End  grid  point  spacing  factor 

0.09 

Surface-normal  stretching  type  at  nose 

Exponential 

Surface-normal  stretching  type  at  end 

Exponential 

Elliptical  smoothing  iterations 

3 

Overlap  correction  factor 

0.30 

Normal  relaxation  factor 

0.30 

Segment 

Angle  (deg) 

Length  (m) 

Num  Points 

Spacing 

Lower  Segment  1 

7 

0.0267 

150 

Automatic 

Lower  Segment  2 

7 

0.53808 

211 

Automatic 

Upper  Segment  1 

35 

0.0003 

10 

Automatic 

Upper  Segment  2 

32.5 

0.0002 

10 

Automatic 

Upper  Segment  3 

30 

0.0005 

10 

Automatic 

Upper  Segment  4 

25 

0.0007 

10 

Automatic 

Upper  Segment  5 

20 

0.0012 

10 

Automatic 

Upper  Segment  6 

18 

0.0025 

10 

Automatic 

Upper  Segment  7 

16 

0.006 

10 

Automatic 

Upper  Segment  8 

14 

0.03 

10 

Automatic 

Upper  Segment  9 

13.1 

0.6 

281 

Automatic 
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D.2  PSE  Script 

# ! /usr/bin/perl 

require  . /lib/pse_scripting . lib" ; 
my  $od=  "/home/trobarge/pseruns/139" ; 
my  $pid=  $$; 
my  $np=2; 

##  Locations  at  which  to  start 

##  Want  to  be  able  to  specify  x  distances  to  match  experiments 

##  and  automatically  convert  them  to  s  distances 

##  Use  general  formula  for  x  to  s  on  frustum  of  sphere-cone 

#  Set  up  parameters : 

my  $Rn=0 . 000508;  #  Nose  radius  in  meters 
my  $phi=7;  #  Cone  half  angle  in  degrees 

my  @xdist=  map  {$_/2;>  (1..16);  #  Distances  to  run  in  inches 
my  (@sdist) ; 

#  Convert  phi  to  radians : 
my  $pi=3. 14159; 

my  $phirad=$phi*$pi/180 ; 

#  Convert  each  point 
foreach  my  $x  (@xdist)  { 

my  $slong=$pi*$Rn/ 180* (90-$phi) + 

($x* . 0254-$Rn* (1-sin ($phirad) /cos ($phirad) ) ) /cos ($phirad) ; 
$slong=~/([\d\.]{6»/; 
push  @sdist,  $1; 

} 

##  Frequencies  to  run 
my  @freqs=  map  {$_*10;}  (1..30); 

#  Set  PSE-Chem  parameters 
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my  °/0params=  (iest_omeg=>0 , iftest3=>l , lstest3=>0 , 
f low=>0 . 2 , f high=>2 . 5 , istep_march=>3 , 
mean_f low_f ile=>" ’ . . /out/shannl50 .bin, ") ; 

#  Make  the  main  output  directory  if  it  doesn't  already  exist 
if  (!-d  $od)  { 
mkdir($od) ; 
if  ( ! -d  $od)  i 

die  "Unable  to  make  output  directory  ’ $od’\n"; 

> 

> 

foreach  my  $f  (@freqs)  { 

$params{"omegal"}=  $f;  $params{"omega2"}=  $f ; 
my  $count=0; 

foreach  my  $s  (@sdist)  { 

#  Identify  runs  by  xdist  in  inches 
my  $run="$f-$xdist [$count] " ; 

$params{"start_sdist"}=  $s; 
my  @oargs=  (-np=>$np,-copy_to=>$od, 

-args=>">  $od/psechem. $run . out" , 
-id=>$run,-outdir=>" . ,/out/pserun.$pid.data" , 

-machinef ile=>$machines) ; 

#  Run  PSE-Chem 

print  " —  Running  PSE-Chem  on  run  $run\n" ; 

my  °/0ret=  psechem_run_solver (-inpf  ile=>\°/0params ,  @oargs)  ; 

psechem_check_error (%ret) ; 

$count++ ; 


} 

exit ; 
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Table  D.2:  Selected  mean-flow  data  points  for  Rufer’s  45  psia  blunt  cone. 


i 

3 

X 

y 

u 

V 

T 

rhol 

100 

100 

0.000649426 

0.00054862 

166.748 

23.6314 

302.1 

0.00851877 

200 

100 

0.00556437 

0.00120189 

169.859 

20.8005 

304.759 

0.00401514 

300 

100 

0.153165 

0.0200131 

475.453 

58.0136 

255.638 

0.00494973 

400 

100 

0.423942 

0.0545917 

761.355 

93.5567 

138.664 

0.00904155 

380 

1 

0.361808 

0.0448739 

0 

0 

300 

0.00418237 

380 

21 

0.361782 

0.0450876 

90.1133 

11.0512 

304.538 

0.00412025 

380 

41 

0.361748 

0.0453598 

203.443 

24.9148 

301.709 

0.00415896 

380 

61 

0.361707 

0.0456987 

343.936 

42.0663 

284.934 

0.00440401 

380 

81 

0.361655 

0.0461207 

517.583 

63.2922 

243.29 

0.00515838 

380 

101 

0.36159 

0.046646 

719.269 

88.2518 

161.753 

0.00776011 

380 

121 

0.36151 

0.0472999 

848.891 

104.766 

72.5221 

0.0173031 

380 

141 

0.36141 

0.0481142 

851.67 

103.487 

67.8951 

0.0184762 

380 

161 

0.361286 

0.0491279 

851.916 

101.358 

67.9034 

0.0184585 

380 

181 

0.361131 

0.0503899 

852.243 

98.8221 

67.8784 

0.0184351 
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Table  D.3:  Selected  PSE  N-factor  data  for  Rufer’s  45  psia  blunt  cone. 


LO 

(kHz) 

Starting  x 

(inches) 

s  =  0.1 

s  =  0.2 

N  factor 

s  =  0.3 

s  =  0.4 

s  =  0.5 

100 

4.0 

0 

0.0485344 

0.1259 

0.240659 

0.546947 

120 

5.0 

0 

0.0226291 

0.159789 

0.920884 

2.69965 

140 

1.0 

0 

0.0422786 

0.846723 

2.76787 

3.2411 

160 

6.5 

0 

0.195131 

2.03906 

2.30514 

1.33935 

180 

0.5 

0 

0.988848 

1.9931 

0.857645 

0 

200 

4.0 

0 

1.50983 

0.74352 

0 

0 

220 

1.0 

0.118642 

1.20406 

0 

0 

0 

240 

1.0 

0.333265 

0.298658 

0 

0 

0 

260 

2.5 

0.493529 

0 

0 

0 

0 

280 

2.0 

0.582963 

0 

0 

0 

0 

300 

1.5 

0.288483 

0 

0 

0 

0 

